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Abstract 


Laser  Radar  sensors  can  be  designed  to  provide  two-dimensional  (2-D)  and  three- 
dimensional  (3-D)  images  of  a  scene  from  a  single  laser  pulse.  Currently,  there  are 
various  data  recording  and  presentation  techniques  being  developed  for  3-D  sensors. 
While  the  technology  is  still  being  proven,  many  applications  are  being  explored 
and  suggested.  As  technological  advancements  are  coupled  with  enhanced  signal 
processing  algorithms,  it  is  possible  that  this  technology  will  present  exciting  new 
military  capabilities  for  sensor  designers  and  end  users. 

The  goal  of  this  work  is  to  develop  an  algorithm  to  enhance  the  utility  of  3-D  Laser 
Radar  sensors  through  accurate  ranging  to  multiple  surfaces  per  image  pixel  while 
minimizing  the  effects  of  diffraction.  Via  a  new  3-D  blind  deconvolution  algorithm, 
it  will  be  possible  to  realize  numerous  enhancements  over  both  traditional  Gaussian 
mixture  modeling  and  single  surface  range  estimation.  While  traditional  Gaussian 
mixture  modeling  can  effectively  model  the  received  pulse,  we  know  that  its  shape 
is  likely  altered  due  to  optical  aberrations  from  the  imaging  system  and  the  medium 
through  which  it  is  imaging.  Simulation  examples  show  that  the  multi-surface  rang¬ 
ing  algorithm  derived  in  this  work  improves  range  estimation  over  standard  Gaussian 
mixture  modeling  and  frame-by-frame  deconvolution  by  up  to  89%  and  85%  respec¬ 
tively. 
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IMPROVING  MULTIPLE  SURFACE  RANGE  ESTIMATION  OF  A 
3-DIMENSIONAL  FLASH  LADAR  IN  THE  PRESENCE  OF  ATMOSPHERIC 

TURBULENCE 

I.  Introduction 

Three  dimensional  FLASH  LAser  Detection  and  Ranging  (LADAR)  sensors  are 
a  special  class  of  Light  Detection  and  Ranging  (LIDAR)  sensors  that  are  able  to 
provide  precise  range  measurements  for  every  pixel  in  an  imaged  scene.  Interest  in 
Three-Dimensional  (3-D)  FLASH  LADAR  systems  is  increasing  for  both  military 
and  civilian  applications  over  the  more  traditional  scanning  LADAR  system.  This 
increase  in  popularity  is  due  to  the  fact  that  3-D  FLASH  LADAR  systems  can  obtain 
an  entire  3-D  image  or  data  cube  with  a  single  laser  pulse.  In  an  ideal  environment, 
the  spatial  resolution  that  this  class  of  sensors  can  achieve  is  limited  only  by  the  ability 
to  design  and  manufacture  components  with  precision.  However,  in  an  operational 
environment,  the  images  may  be  distorted  not  only  by  system  limitations,  but  also 
by  the  atmosphere  through  which  the  light  must  pass. 

1.1  Motivation 

The  motivation  for  this  research  initially  stemmed  from  an  inquiry  from  the  De¬ 
partment  of  Homeland  Security  (DHS)  about  imaging  through  obscurations.  Current 
state-of-the-art  tactical  sensors  such  as  the  AN/AAQ  33  -  SNIPER  targeting  pod 
shown  in  Figure  1.1,  rely  on  passive  target  illumination.  While  this  type  of  illumi¬ 
nation  has  many  inherent  advantages,  it  does  introduce  several  limitations  as  well. 
Like  all  passive  Electro  Optical  /  Infra-Red  (EO/IR)  sensors,  SNIPER  suffers  from 
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Figure  1.1:  AN/AAQ  33  -  SNIPER  targeting  pod  integrated  onto  the  B-1B. 

thermal-crossover  associated  with  the  natural  heating  and  cooling  of  the  earth’s  sur¬ 
face.  Further,  EO/IR  sensors  rely  on  the  addition  of  an  active  illumination  designator 
to  perform  ranging  to  a  target.  In  addition  to  tactical  sensors,  a  study  produced  by 
Forecast  International  in  2011  estimated  that  nearly  $17  billion  dollars  will  be  spent 
on  new  remote  sensing  satellite  technology  between  2012  and  2021  [1],  LADAR  tech¬ 
nology  is  one  of  the  relatively  new  imaging  technologies  that  will  likely  be  employed 
on  the  next  generation  of  remote  sensing  satellite  platforms. 

The  aforementioned  are  just  a  sampling  of  the  justification  for  employing  3-D 
LADAR  technology  on  future  generation  remote  sensors.  Unfortunately,  due  to  con¬ 
straints  such  as  low  spatial  resolution  and  limits  on  laser  power,  the  current  state  of 
FLASH  LADAR  technology  would  yield  limited,  if  any  improvement  in  total  capabil¬ 
ity  over  the  currently  fielded  passive  sensors.  However,  it  may  be  possible  to  develop 
algorithms  that  limit  the  negative  impacts  of  the  operational  environment  and  couple 
them  with  technological  advances  in  sensor  resolution  and  sensitivity  to  yield  vast 
improvements  in  future  imaging/targeting  sensors. 

Due  to  the  employment  of  active  illumination,  3-D  FLASH  LADAR  sensors  can 
gather  ranging  information  for  every  point  within  a  targeted  scene  nearly  simulta- 
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neously.  Further,  depending  on  the  detection  methodology  employed,  the  possibility 
of  imaging  through  obscurations  becomes  a  reality.  Imaging  through  obscurations 
has  numerous  civilian  and  military  applications.  For  instance,  Advanced  Scientific 
Concepts  (ASC),  the  manufacturer  of  the  sensor  that  will  be  utilized  throughout  this 
research  effort  is  currently  interested  in  demonstrating  the  ability  to  detect  targets 
obscured  by  smoke,  fog  or  dust.  One  particular  investigation  underway  applies  to  the 
tracking  of  a  refueling  drogue  by  an  autonomous  aerial  vehicle  [13].  The  images  in 
Figure  1.2  demonstrate  a  valuable  capability  inherent  with  3-D  FLASH  LADAR. 


(a)  Original  Image  (b)  Target  Visually  Obscured  by  Smoke 


(c)  First  Return  Detection  (d)  Second  Return  Detection 


Figure  1.2:  In  this  figure  we  show  a  scene  visually  obscured  by  smoke.  We  then  show 
that  the  3-D  FLASH  LADAR  system  is  able  to  see  through  the  smoke  and  detect  the 
target  of  interest. 
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Here  we  demonstrate  that  the  LADAR  system  is  easily  able  to  detect  targets 
that  are  visually  obstructed  by  smoke.  Figures  1.2(a)  and  1.2(b)  were  taken  with  a 
traditional  video  camera  while  Figures  1.2(c)  and  1.2(d)  were  taken  with  a  3-D  FLASH 
LADAR  system.  The  3-D  FLASH  LADAR  system  does  receive  a  return  off  of  the 
smoke  as  shown  in  Figure  1.2(c).  However,  the  multi-surface  ranging  capability  that 
this  dissertation  will  focus  on  allows  us  to  detect  additional  reflections,  thus  revealing 
the  targets  of  interest  in  Figure  1.2(d).  Clearly  this  has  numerous  military  as  well  as 
civilian  applications.  On  the  civilian  side,  this  technology  has  applications  for  vehicles 
traveling  through  fog  or  for  firefighters  trying  to  see  through  a  smoky  room  to  rescue 
trapped  personnel.  On  the  military  side,  this  technology  would  be  useful  for  areas 
such  as  landing  an  aircraft  or  helicopter  in  brown  out  conditions  or  conducting  aerial 
refueling  operations  through  clouds  among  others. 

Multi-surface  ranging  with  the  use  of  3-D  FLASH  LADAR  can  also  be  useful  in 
accurately  discriminating  camouflaged  targets  of  interest.  A  tactic  that  is  commonly 
employed  on  the  battlefield  is  to  develop  mock  targets  of  low  value.  These  targets 
are  difficult  to  discern  in  a  tactical  environment  where  the  attack  may  be  conducted 
from  miles  away  in  a  fast  moving  vehicle.  Thus,  it  can  be  an  effective  tactic  because 
it  forces  the  aggressor  to  use  a  potentially  high  valued  weapon  on  a  target  of  little 
to  no  importance.  In  the  visual  spectrum  it  is  often  easy  to  generate  false  targets  or 
hide  actual  targets  with  camouflage  netting.  However,  the  ability  of  a  3-D  FLASH 
LADAR  system  to  detect  images  with  range  diversity  can  make  it  far  more  difficult 
to  design  false  targets  or  effectively  camouflage  real  ones. 

Clearly,  3-D  FLASH  LADAR  provides  additional  useful  information  to  the  user. 
In  a  traditional  2-D  image,  we  only  detect  some  sort  of  intensity  information  in  the 
spatial  domain.  For  3-D  FLASH  LADAR  images,  we  not  only  receive  intensity  and 
range  information,  but  depending  on  the  target  geometry  or  physical  characteristics, 
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we  may  also  be  able  to  discern  information  from  the  width  of  the  reflected  pnlse. 
This  pulse  width  information  can  often  be  used  to  discern  whether  a  target  is  diffuse 
or  solid.  Additionally,  it  may  be  possible  to  obtain  target  orientation  based  on  pulse 
width  expansion  [26]. 

Ultimately,  the  ability  to  accurately  assess  the  threat  environment  is  critically 
important  for  numerous  reasons.  Detecting  targets  that  may  be  concealed  by  man¬ 
made  camouflage  or  natural  obscuration  has  long  been  a  goal  of  numerous  sensor 
developments,  as  those  targets  may  pose  a  threat  to  the  safety  of  forces  or  success 
of  a  military  operation.  Optical  sensors  based  on  active  laser  illumination  such  as 
FLASH  LADAR  have  the  distinct  advantage  of  operating  at  extremely  short  wave¬ 
lengths  enabling  light  to  pass  through  small  voids  in  the  obscuration  and  reflect  off  of 
potential  surfaces  of  interest.  The  ability  to  finely  sample  the  returned  waveform  with 
a  high  resolution  detector  array  will  enable  an  imaging  sensor  to  produce  an  accurate 
representation  of  the  obscured  target  [50].  In  addition  to  the  previously  mentioned 
applications,  advancement  of  this  technology  has  potential  application  in  areas  such 
as  terrain  mapping,  forestry  classification  and  autonomous  vehicle  navigation. 

Current  manufacturing  limitations  for  3-D  FLASH  LADAR  are  often  centered 
around  sampling  rates,  detector  size  and  development  of  optical  components  that 
are  free  from  aberration.  The  dimensions  of  each  pixel  in  the  detector  array  are 
currently  the  primary  limiting  factor  on  attainable  spatial  resolution.  As  with  all 
emerging  camera  technologies,  we  expect  this  to  improve  as  manufacturing  processes 
are  refined.  Unfortunately,  as  a  consequence  of  nonuniform  heating  and  cooling  of 
the  Earth’s  atmosphere,  the  temperature- induced  inhomogeneities  of  the  refractive 
index  of  the  air  may  have  a  significant  impact  on  attainable  spatial  resolutions  [21]. 
This  dissertation  will  address  those  problems  and  reinforce  that  previously  ill-posed 
problems  for  2-D  imaging  may  actually  be  overdetermined  for  3-D  imaging  [49]. 
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1.2  LADAR  Technology 
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Figure  1.3:  LADAR  flow  chart  with  the  specific  technology  that  the  bulk  of  this  research 
will  be  concerned  with  down  the  right  hand  side  highlighted  in  green. 

LAser  Detection  and  Ranging  (LADAR)  systems  are  a  subset  of  Light  Detection 
and  Ranging  (LIDAR)  systems  where  the  delineation  is  tied  to  the  type  of  illumination 
source.  Further  exploration  of  this  class  of  sensor  warrants  a  brief  description  of  the 
numerous  variations  on  the  technology.  Figure  1.3  encapsulates  many  of  the  numerous 
variants  of  LADAR  technology.  There  are  currently  two  widely  recognized  forms  of 
3-D  LADAR  systems.  One  type  is  a  scanning  LADAR  system  where  each  pixel  in  the 
image  is  a  result  of  measuring  the  return  from  a  separate  laser  pulse.  The  scanning 
variant  of  LADAR  has  been  widely  studied,  and  many  algorithms  have  been  developed 
to  process  the  measured  data  such  as  in  [31],  [32]  and  [63].  While  the  Airborne  Laser 
Scanner  ALS  variants  of  3-D  LADAR  have  demonstrated  unique  capability  in  remote 
sensing,  they  commonly  require  a  significant  amount  of  time  to  form  a  complete 


6 


image.  The  time  delay  associated  with  this  technique  will  likely  be  unacceptable  for 
employment  on  a  fast  moving  platform  operating  in  a  highly  dynamic  environment. 
Fortunately,  continuing  technical  advancements  are  giving  rise  to  systems  known  as 
3-D  FLASH  LADAR  systems. 

Unlike  scanning  LADAR  systems,  a  3-D  FLASH  LADAR  is  able  to  form  a  com¬ 
plete  3-D  image  or  data  cube  by  simultaneously  measuring  the  returned  pulse  for 
every  pixel  in  an  image.  The  3-D  image  or  data  cube  is  essentially  a  series  of  images 
where  an  Avalanche  Photo  Diode  (APD)  array  measures  the  returning  photons  for 
each  pixel  separated  by  a  constant  time  interval.  Figures  1.4  and  1.5  provide  a  sim¬ 
plified  sketch  of  a  3-D  FLASH  LADAR  in  operation.  Commonly,  a  laser  illuminator 


3D  FLASH  LADAR 


Figure  1.4:  Transmit  Portion  of  3-D  FLASH  LADAR  Operation. 

3D  FLASH  LADAR 


Figure  1.5:  Receive  Portion  of  3-D  FLASH  LADAR  Operation. 
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will  fire  a  short  pulse  of  light  through  some  sort  of  beam  spreader  to  achieve  the  de¬ 
sired  illumination  of  the  target  area.  In  a  full-waveform  system,  the  APD  array  will 
then  simultaneously  sample  the  reflected  pulse  for  each  pixel  in  the  detector  array 
at  a  predetermined  rate.  A  separate  2-D  image  will  be  formed  for  each  sampling  of 
the  reflected  laser  pnlse.  Based  on  this  methodology,  accurate  ranging  will  be  highly 
dependent  upon  the  ability  to  precisely  account  for  the  time  of  flight  of  the  laser 
pulse. 

The  sampling  rate  and  number  of  samples  within  the  range  gate  will  be  largely 
system/mission  dependent;  however,  it  is  worthwhile  to  consider  the  two  common 
methods  for  triggering  the  start  and  end  of  the  range-gate.  The  first  mode  under 
consideration,  Staring  Underwater  LAser  Radar  (SULAR)  mode,  records  the  first 
frame  of  data  at  a  pre-dehned  time  from  when  the  pulse  was  fired.  In  this  manner, 
the  range  gate  for  a  3-D  image  will  have  a  fixed  window  for  each  detector  based  on 
the  start  time,  time  between  frames  and  total  number  of  frames  in  the  3-D  image. 
Another  mode  of  operation  is  commonly  referred  to  as  “HIT  mode”  where  each  image 
pixel  may  trigger  the  start  of  the  range  gate  at  a  different  time  based  on  the  received 
signal  level.  The  bulk  of  the  research  presented  in  this  dissertation  will  focus  on  3-D 
FLASH  LADAR  operating  in  SULAR  mode. 

1.3  Research  Contributions 

The  following  subsections  offer  a  brief  description  of  each  of  the  three  core  areas 
of  research  covered  in  this  dissertation  as  well  as  their  associated  contributions. 


1.3.1  Parameterized  Blind  Deconvolution  Through  Convergence  of 
Variance  (Chapter  III). 

This  area  of  research  first  revisits  algorithms  originally  presented  by  MacManus 
and  MacDonald  for  accomplishing  parameterized  blind  deconvolution  of  Two  -  Di¬ 
mensional  (2-D)  images  in  the  presence  of  Poisson  and  negative  binomial  noise  [41], 
[45].  The  Convergence  of  Variance  (CoV)  algorithm  developed  by  MacManus  and 
the  Maximum  a-priori  (MAP)  estimate  developed  by  MacDonald  each  offer  blind 
deconvolution  methods  tailored  for  implementation  in  a  tactical  environment  where 
processing  time  is  extremely  important.  Initially,  a  comparison  of  the  two  algorithms 
is  offered  and  several  unique  conclusions  are  developed.  Ultimately  a  new  explanation 
is  given  which  reveals  that  the  MAP  estimate  provides  no  new  capability  over  the 
CoV  technique.  Further,  an  extension  to  the  CoV  algorithm  is  proposed  in  order 
to  perform  parameterized  blind  deconvolution  on  3-D  FLASH  LADAR  images.  The 
primary  goal  of  this  portion  of  research  was  to  demonstrate  the  ability  to  find  a  pa¬ 
rameterized  Point  Spread  Function  (PSF)  which  could  be  used  in  the  multi-surface 
ranging  algorithm  presented  in  Chapter  IV. 

1.3.2  Multiple  Surface  Detection  and  Estimation  (Chapter  IV). 

The  core  area  of  research  for  this  dissertation  was  the  development  of  a  novel  itera¬ 
tive  algorithm  using  an  Expectation  Maximization  (EM)  strategy  which  will  simulta¬ 
neously  solve  for  multiple  ranges  per  image  pixel  and  remove  the  effects  of  diffraction 
[14],  Challenges  with  employing  an  EM  strategy  are  discussed,  and  the  techniques 
developed  to  mitigate  those  challenges  are  presented.  Additionally,  a  constraint  is  pro¬ 
posed  and  applied  to  the  initially  developed  multi-surface  ranging  algorithm  which 
further  enhances  its  capability.  Both  the  constrained  and  non-constrained  multi¬ 
surface  algorithms  which  account  for  the  effects  of  diffraction  are  then  compared 
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against  an  algorithm  which  does  not  account  for  diffraction  and  one  that  employs 
traditional  2-D  image  deconvolution  techniques  to  account  for  diffraction.  This  new 
algorithm  is  ultimately  shown  to  provide  a  significant  improvement  in  multi-surface 
ranging  capability.  This  capability  will  be  critical  in  employing  3-D  FLASH  LADAR 
technology  in  an  environment  where  the  ability  to  image  through  obscurations  is 
required. 

1.3.3  Imaging  Through  Turbulence  (Chapter  V). 

Finally,  the  next  step  in  enhancing  the  multi-surface  detection  algorithm  presented 
in  Chapter  IV  is  developed.  Initially  the  multi-surface  ranging  algorithm  assumed  a 
known  PSF  to  account  for  the  effects  of  diffraction.  Determining  the  PSF  is  a  chal¬ 
lenging  problem  known  as  blind  deconvolution  that  has  been  the  focus  of  a  wealth 
of  research.  Traditionally,  blind  deconvolution  problems  for  2-D  images  are  ill-posed. 
However,  this  research  shows  that  a  maximum  likelihood  approach  to  estimate  the 
PSF  parameterized  by  Fried’s  seeing  parameter,  r0,  is  possible  with  the  addition  of 
range  diversity  from  3-D  images.  Unlike  the  MAP  estimate  proposed  by  MacDonald 
and  discussed  in  Chapter  III,  this  technique  does  not  require  the  introduction  of  a 
prior  distribution  for  the  value  of  r0.  Alternatively,  the  CoV  technique  developed  by 
MacManus  could  first  be  used  to  identify  a  parameterized  PSF  before  proceeding  to 
use  the  multi-surface  ranging  algorithm  to  develop  the  surface  profile.  However,  this 
research  will  show  that  both  the  value  for  ro  and  the  surface  profile  can  be  estimated 
simultaneously.  Systems  of  equations  will  be  provided  that  highlight  the  added  capa¬ 
bility  of  working  this  type  of  problem  with  3-D  FLASH  LADAR  images.  Finally,  the 
results  will  be  verified  through  the  use  of  both  simulation  and  experimental  data. 
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1.4  Organization 


This  document  will  be  organized  as  follows.  Relevant  background  material  on 
LADAR  technology,  theoretical  models  and  experimental  data  collection  systems  will 
be  presented  in  Chapter  II.  Chapter  III  examines  two  previously  developed  parameter¬ 
ized  blind  deconvolution  algorithms,  and  discusses  their  application  to  3-D  FLASH 
LADAR  imagery.  Chapter  IV  will  detail  the  derivation  of  a  multi-surface  ranging 
algorithm  that  simultaneously  develops  the  range  profile  and  removes  the  effects 
of  diffraction  from  a  3-D  FLASH  LADAR  image,  given  a  known  PSF.  Chapter  V 
presents  a  joint  estimation  technique  for  simultaneously  developing  the  range  profile 
and  parameterized  PSF.  Finally,  Chapter  VI  will  summarize  the  conclusions  from  this 
research  and  present  ideas  for  future  research. 
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II.  Background 


This  chapter  presents  necessary  support  material  for  the  research  presented  in  this 
dissertation.  First,  a  model  will  be  presented  for  simulation  and  representation  of  the 
data  received  by  a  3-D  FLASH  LADAR  sensor.  Additionally,  theory  will  be  presented 
to  identify  the  ill-effects  of  imaging  through  a  turbid  medium  in  Section  2.1.  Common 
techniques  for  improving  the  overall  Signal  to  Noise  Ratio  (SNR)  in  the  detected 
images  will  then  be  presented  in  Section  2.2.  Section  2.3  will  focus  on  a  2-D  and  3-D 
sensor  that  will  be  used  throughout  the  research  effort  for  the  collection  of  measured 
data.  An  important  concept  for  any  iterative  algorithm  is  how  to  determine  when  the 
algorithm  can  terminate.  This  research  will  employ  the  convergence  of  variance  as 
the  stopping  criteria  as  discussed  in  Section  2.4.  In  simulation,  the  blur  applied  to  the 
detected  image,  or  PSF,  is  known.  However,  with  experimental  data,  we  must  be  able 
to  measure  the  PSF  for  comparison  with  the  estimates  obtained  from  the  algorithms 
developed.  Section  2.5  will  present  the  method  used  for  measuring  Fried’s  seeing 
parameter,  ro,  which  will  be  used  in  the  parameterized  model  for  the  PSF.  Finally, 
this  chapter  will  present  a  summary  of  similar  research  previously  accomplished  in 
Section  2.6. 

2.1  3-D  FLASH  LADAR  Theory 

Three  dimensional  FLASH  LADAR  is  fundamentally  different  from  scanning  vari¬ 
ants  in  that  the  entire  remote  scene  is  illuminated  with  a  single  short  pulse  of  the  laser 
illuminator.  On  the  other  hand,  the  scanning  variants  only  illuminate  a  narrow  Field 
of  View  (FOV)  corresponding  with  a  single  image  pixel  per  pulse  of  the  laser.  This 
narrow  beam  is  then  scanned  over  the  target  area  to  compile  an  entire  3-D  image. 
While  scanning  systems  have  been  shown  to  provide  extremely  detailed  images,  the 
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size,  complexity,  low  frame  rate  and  cost  of  stabilization  and  beam  steering  hardware 
make  the  technology  impractical  for  tactical  use. 

This  research  will  utilize  some  concepts  learned  through  previous  research  with 
scanning  variants  of  LADAR;  however,  the  research  will  focus  on  the  3-D  FLASH 
LADAR  variants.  Since  the  entire  remote  scene  is  illuminated  with  a  single  pulse 
of  the  laser,  the  resultant  data  collected  will  depend  on  the  triggering  mechanism 
employed.  As  shown  in  Figure  1.3,  the  data  stored  is  primarily  classified  as  either 
Range  and  Amplitude  (R&A)  or  full  waveform.  For  the  R&A  data,  each  pixel  in  the 
detector  will  employ  some  sort  of  detection  methodology  to  determine  the  range  and 
intensity  from  the  reflected  waveform.  The  only  data  stored  for  each  pixel  will  be 
the  resultant  range  and  amplitude.  In  full  waveform  data,  there  are  various  forms 
of  triggering  mechanisms  worthy  of  discussion.  The  ASC  systems  employ  either  HIT 
mode  or  SULAR  mode  of  triggering.  HIT  mode  allows  for  each  pixel  within  the 
detector  to  have  a  unique  range  gate  where  the  start  of  the  range  gate  is  triggered  at 
some  intensity  threshold.  While  the  HIT  mode  of  operation  definitely  has  application, 
the  techniques  presented  in  this  research  require  that  each  pixel  in  the  image  have 
the  same  range  gate. 

For  the  SULAR  mode  of  operation,  each  pixel  within  the  detector  will  have  pre¬ 
cisely  the  same  range  gate.  Each  time  the  detector  is  sampled,  a  separate  2-D  frame 
will  be  stored  which  contains  an  intensity  representation  corresponding  with  the  num¬ 
ber  of  photons  received  over  the  detector  integration  time.  Therefore,  each  2-D  frame 
will  also  have  a  discrete  range  corresponding  to  its  sample  time.  This  compilation 
of  data  frames  can  then  be  stacked  into  a  3-D  data  cube  such  that  each  pixel  in  the 
detector  will  have  the  reflected  waveform  corresponding  to  the  target  area  and  the 
established  range  gate.  Based  on  this  mode  of  operation,  it  is  possible  that  certain 
portions  of  the  scene  may  not  result  in  a  pulse  return  within  the  range  gate.  These 
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non-illuminated  areas  of  the  scene  will  therefore  only  measure  the  detector  bias  for 
each  sample  within  the  range  gate. 

Accounting  for  range  once  a  certain  threshold  is  achieved  or  where  a  peak  is 
detected  can  only  provide  range  accuracy  to  the  nearest  sample.  Provided  the  sample 
rate  is  high  enough,  full  waveform  data  storage  allows  for  highly  precise  ranging 
capability.  Generally  we  would  want  to  sample  the  detector  at  the  Nyquist  rate. 
Then,  using  techniques  such  as  the  correlation  method  presented  by  Richmond  and 
Cain  [58],  or  the  EM  algorithm  presented  in  this  research,  we  can  achieve  sub-sample 
ranging  accuracy. 


2.1.1  Pulse  Model. 

This  research  will  employ  the  common  technique  of  modeling  the  transmitted  and 
received  pulse  for  a  LADAR  system  as  a  Gaussian  (single-surface  case)  or  mixture  of 
Gaussians  (multi-surface  case).  While  the  outgoing  waveform  will  largely  dictate  the 
inaccuracy  induced  by  the  Gaussian  approximation,  for  many  sensors  this  model  is 
generally  considered  acceptable  [46] .  Future  work  may  consider  the  error  introduced; 
however,  the  error  is  likely  to  be  sensor  dependent  [58].  Using  an  early  variant  of  the 
FLASH  3-D  LADAR  sensor  developed  by  ASC,  the  Gaussian  approximation  allowed 
for  accurate  single  surface  ranging  algorithms  to  be  developed  by  Dolce  [16]  and 
McMahon  [48] .  The  Gaussian  model  for  the  received  intensity  of  the  pulse  P  (t)  is 


p(t) 


e  l 

2  a2 


(2.1) 


where  t  is  time,  A  is  the  pulse  amplitude  and  a  is  the  width  parameter. 

Using  the  Gaussian  model  as  our  basis,  much  of  this  dissertation  will  consider  the 
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possibility  of  N  surface  returns  per  image  pixel  according  to 


N 

P  (x,  y,  rk)  =  P[n)  {x,  V,  rk).  (2.2) 

n= 1 


Where  each  of  the  N  individual  Gaussian  pulses,  P^n\  reflected  by  the  target  at  a 
range  of  r,(n)  is 


P(n)  (x,  y,  rk) 


A ^  ( x ,  y)  \  (rk  ~  r(n)  ( x ,  y))2 

, _  ,  , - exp - p- 

( x ,  y)  2  (cdn)  (x,  y)) 


(2.3) 


The  indexes  x  and  y  indicate  the  area  in  the  target  plane  corresponding  to  the  indi¬ 
vidual  pixels  in  the  M  x  M  detector  array,  and  the  index  rk  represents  the  discrete 
range  for  the  kth  frame  in  the  data  cube,  or  in  other  words,  each  time  the  pulse  is 
sampled  by  the  sensor.  The  nth  amplitude,  pulse  width  and  range  of  the  pulse  mixture 
are  represented  by  A^n\  al'n'>  and  r ^  respectively. 

The  received  pulse  is  modeled  by  the  intensity  function  I(u,v,  rk),  where  the 
received  intensity  is  found  by  convolving  the  pulse  and  the  PSF,  h,  as  shown  in 


M  M 

I  (■ u,v,rk )  =  EE  P  (x,  y,  rk)  h(u  —  x,v  —  y),  (2.4) 

x=l  y= 1 


where  the  indexes  u  and  v  represent  the  detector  plane  coordinates.  It  has  been 
demonstrated  that  the  number  of  photons  that  arrive  during  the  detector’s  integration 
time  can  be  modeled  with  Poisson  statistics  [21].  Section  2.2  will  discuss  the  choice 
of  Poisson  statistics  over  negative  binomial  statistics  in  more  detail.  In  addition  to 
the  laser  light  reflected  off  of  the  target,  the  detector  may  also  receive  photons  from 
background  lighting,  and  thermal  noise.  These  additional  photon  sources  will  be 
accounted  for  as  a  pulse  bias  for  each  detector  in  the  APD  array  B  (u,v).  The  noise 
generated  by  the  bias  will  also  follow  a  Poisson  distribution,  resulting  in  a  model  for 
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the  total  intensity  measured  by  the  detector  of 


hot  (■ u,v,rk ) 


'  N  M  M 

EEE  P(re)  (x,  y,  rk)  h  ( u 
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x,  v  -  y) 


+  B  (u,v) . 


(2.5) 


In  Section  4.1,  a  method  is  derived  for  estimating  the  pulse  bias  in  conjunction  with 
the  amplitude,  pulse  width  and  range  of  the  pulse.  In  reality,  the  signal  bias  could  be 
measured  by  producing  an  image  in  the  absence  of  a  laser  pulse;  however,  this  may 
not  always  be  possible,  thus  driving  the  additional  need  to  estimate  this  parameter. 

If  we  assume  independence  of  the  measurements  at  each  time  step  and  for  every 
pixel  in  the  detector  array,  we  can  state  the  joint  probability  of  the  observed  data,  d, 
as 
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(2.6) 


To  minimize  confusion  in  notation,  an  upper-case  P  will  be  used  to  represent  the  pulse, 
and  lower-case  p  will  be  used  in  representation  of  the  various  probabilities  throughout 
this  work.  The  expected  value  of  the  noisy  2-D  frame,  Dk(u,v),  corresponding  with 
a  range  to  target  of  rk  is  hot(u , v ,  fk). 

In  summary,  each  2-D  frame  is  impacted  by  the  effects  of  diffraction,  an  additive 
bias  and  noise  as  shown  in  Figure  2.1.  The  noisy  received  3-D  image,  d,  of  a  remote 
object,  o,  is  comprised  of  a  series  of  2-D  frames.  The  received  intensity  for  each  frame, 
hot,  is  the  summation  of  an  additive  bias  and  the  blurred  image,  i. 
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d 


Figure  2.1:  Model  for  received  data  that  is  impacted  by  diffraction  and 


noise. 


17 


2.1.2  Spatial  Sampling  Requirements. 


In  the  diffraction  limited  case,  the  maximum  achievable  spatial  frequency,  z/maa;, 
can  be  computed  by 


D 

hnax  ,  j,  > 

AJl 


(2.7) 


where  A  is  the  wavelength  of  the  light  of  interest  and  fi  is  the  focal  length  of  the  lens 
and  D  is  the  aperture  diameter  [22],  According  to  the  Nyquist  criterion,  the  ability  to 
perfectly  reconstruct  a  digitally  sampled  image,  or  properly  sample,  would  therefore 
require  a  sampling  frequency  of  2//max.  The  importance  of  properly  sampling  an  image 
stems  from  the  use  of  convolution  in  the  model  for  intensity  as  shown  in  (2.5). 

Achieving  Nyquist  sampling  is  a  significant  challenge  in  many  optical  sensing 
applications,  and  especially  in  the  case  of  3-D  FLASH  LADAR.  The  complexity 
of  the  electronics  coupled  with  conventional  optics  often  results  in  a  detector  pixel 
pitch  that  is  much  greater  than  what  is  required  for  proper  sampling.  However, 
assuming  the  standard  progression  of  technology,  history  has  shown  that  the  pixel 
pitch  will  decrease.  The  techniques  presented  in  this  dissertation  are  again  intended 
to  demonstrate  the  potential  value  of  the  methodology.  Future  research  efforts  could 
be  tied  to  exploring  their  utility  where  proper  sampling  may  not  be  possible. 


2.1.3  Sources  of  Image  Blurring  and  Degradation. 

Due  to  numerous  physical  phenomena,  the  images  captured  by  a  sensor  have 
imperfections.  First,  we  know  that  optical  imperfections  with  the  system  can  cause 
the  diffraction  of  light  to  neighboring  areas.  When  operating  a  sensor  within  the 
atmosphere,  optical  imperfections  can  commonly  be  binned  into  those  directly  related 
to  the  manufacturing  of  the  sensor  and  those  that  result  from  atmospheric  turbulence. 
Second,  the  received  image  is  generally  further  degraded  by  the  effects  of  noise.  This 
noise  can  also  result  from  numerous  sources  such  as  read  out  noise,  thermal  noise  and 
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noise  associated  with  the  illumination  source. 


The  total  PSF  or  spatial  impulse  response  of  an  optical  sensor  accounts  for  the 
diffraction  effects  directly  attributed  to  the  sensor  optics,  hoptl  and  those  that  can  be 
attributed  to  atmospheric  turbulence,  hatm-  This  total  PSF,  htot ,  is  the  2-D  convolu¬ 
tion  of  the  two  primary  components  as  shown  in  (2.8),  where  x  and  y  are  the  spatial 
coordinates  of  the  individual  PSFs. 

htot  {x,  y)  =  hopt  (x,  y)  <8>  hatm  (x,  y)  (2.8) 

A  direct  solution  for  the  impulse  response  of  a  non-ideal  aperture  is  difficult  to 
compute;  however,  it  can  be  found  by  conducting  a  propagation  experiment  using 
known  sensor  parameters  [58].  In  frequency  space,  the  total  Optical  Transfer  Func¬ 
tion  (OTF),  Htot,  the  optics  OTF,  Hopt ,  and  atmospheric  OTF,  Hatm,  are  simply  the 
Fourier  transforms  of  their  respective  PSFs,  and  the  convolution  operator  is  replaced 
by  the  multiplication  operator 

{ htot  (Xj  y  ) }  Htot  (pr  j  0/ )  Hopt  (ry ,  z/y)  Hatm  (pr  >  0/ )  >  (2.9) 

where  the  spatial  frequencies  in  two  dimensions  are  parameterized  by  (ux,  uy)  and  T 
is  the  Fourier  operator. 

Imaging  devices  have  an  upper  bound  on  their  cutoff  frequency  dictated  by  their 
optical  specifications  according  to  (2.7).  Imaging  through  a  turbid  medium  will  at 
best  cause  no  further  attenuation  to  the  frequency  content  of  an  image,  but  can  often 
result  in  significant  attenuation.  In  a  2-D  image,  this  attenuation  of  frequency  content 
manifests  as  a  spatial  blurring  of  the  image.  However,  since  the  3-D  FLASH  LADAR 
image  is  a  composition  of  numerous  2-D  frames,  this  spatial  blurring  will  also  cause 
a  temporal  distortion  as  well. 
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As  a  general  rule  of  thumb,  when  collecting  images  with  exposure  times  of  less 
than  jA  of  a  second  we  can  assume  that  the  atmosphere  through  which  the  remote 
scene  is  viewed  remains  constant  [21].  While  the  instantaneous  atmosphere  is  difficult 
to  measure  and  even  more  challenging  to  estimate,  models  exist  to  predict  its  behavior 
on  average.  In  the  short  exposure  scenario  where  the  total  exposure  time  is  less  than 
A  of  a  second,  we  have  the  average  OTF,  HSe, 


a  2/I2P»  +  d)V/'' 

D2  I 

(2.10) 

In  the  mathematical  model  for  Hse  as  a  function  of  spatial  frequency,  the  mean 
wavelength  of  light  detected  is  A  and  r0  is  Fried’s  seeing  parameter.  Applying  this 
model  essentially  has  two  stipulations  of  concern.  First,  we  must  have  a  sufficiently 
short  exposure  time  to  satisfy  the  static  atmosphere  requirement.  Clearly  a  single 
image  taken  with  a  3-D  FLASH  LADAR  system  falls  within  this  requirement.  Second, 
we  must  be  able  to  register  the  images  to  remove  the  motion  effects  of  tip  /  tilt  caused 
by  the  atmosphere.  For  purposes  of  this  research,  we  will  substitute  HSe  as  our  model 
for  the  atmospheric  OTF,  Hatm.  However,  the  techniques  developed  in  subsequent 
chapters  could  also  be  demonstrated  to  work  with  the  simpler  long-exposure  case 
where  the  image  frames  are  averaged  without  motion  compensation.  The  average 
long-exposure  OTF,  Hle,  is 


HSe  (Cr,  vy)  =  exp  <  -3.44 


A  2/,2te  +  *2)' 


5/6 


H 


LE 


\VX  ,Vy, 


=  exp  {  -3.44 


A  2/i2M  +  *S)' 


5/6' 


(2.H) 


Even  though  the  individual  exposures  of  the  3-D  FLASH  LADAR  sensor  are  less  than 
A  of  a  second,  the  use  of  this  model  would  apply  in  cases  where  we  chose  not  to 
account  for  the  tip  /  tilt  motion  caused  by  the  atmosphere. 
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2.2  Image  Variance  and  Effects  of  Averaging 


Due  to  the  partially  coherent  nature  of  the  illumination  source,  the  reflected  pho¬ 
ton  distribution  is  negative  binomial.  This  results  in  a  somewhat  complex  noise 
distribution  for  the  composite  pixel  intensity  for  each  frame.  The  composite  intensity 
also  has  noise  contributions  from  background  illumination  which  is  commonly  mod¬ 
eled  with  a  Poisson  distribution.  Experience  has  shown  that  the  negative  binomial 
noise  from  the  reflected  laser  pulse  is  commonly  the  most  dominant  form  of  noise. 
Therefore  in  regions  of  relatively  high  illumination,  the  other  sources  of  noise  are 
indistinguishable.  The  negative  binomial  Probability  Mass  Function  (PMF)  model 
defines  the  probability  that  we  will  receive  K  photons  over  the  integration  time  as 


p(K) 


r  (k  +  M) 
T(K  +  1)T{M) 


1+^' 

K 


-K 


'  K 
1  +  M 


-M 


(2.12) 


where  M  is  the  coherence  parameter  of  the  laser  illumination  and  K  is  the  expected 
number  of  photons.  Due  to  its  complexity,  the  negative  binomial  PMF  does  not 
lend  itself  to  the  algorithms  that  were  developed  through  the  course  of  this  research. 
Fortunately,  the  Poisson  PMF  has  been  shown  to  be  an  adequate  substitution  based 
primarily  on  the  similarity  between  the  shapes  of  the  two  PMFs  given  the  standard 
range  of  coherence  parameters  [41],  [47].  As  A4  approaches  infinity,  the  variance  of  a 
negative  binomial  random  variable  collapses  to  the  mean, 


a 


2  _ 
nb  ~ 


K 


(2.13) 


thus  converging  to  a  Poisson  variance.  A  comparison  of  the  negative  binomial  and 
Poisson  PMF  is  shown  in  Figure  2.2.  Here  we  see  that  as  JA  increases,  the  negative 
binomial  PMF  does  in  fact  approach  the  Poisson  PMF. 

Motion  Compensated  Frame  Average  (MCFA)  images  are  commonly  used  in  imag- 
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Figure  2.2:  Convergence  of  the  negative  binomial  PMF  to  the  Poisson  PMF  with 
increasing  At  and  K  =  2000. 

ing  applications  to  improve  the  SNR.  The  value  of  using  an  MCFA  can  best  be  un¬ 
derstood  through  a  simple  example.  For  digital  image  processing,  SNR  is  commonly 
modeled  as 

SNR  =  — ,  (2.14) 

a 

or  the  ratio  of  mean  photons  measured  to  the  standard  deviation  of  the  measurement 
[20].  This  definition  of  SNR  is  only  justified  when  the  mean  is  always  non- negative 
such  as  the  case  with  photon  counts.  For  simplicity,  we  will  now  consider  a  Poisson 
process  for  which  we  know  the  mean  is  equal  to  the  variance.  If  our  MCFA  consists  of 
the  summation  of  two  independent  measurements,  both  the  mean  and  variance  will 
increase  by  a  factor  of  two.  However,  the  standard  deviation  only  increases  by  a  factor 
of  y/2.  Therefore,  the  SNR  also  improves  by  a  factor  of  \/2  every  time  we  double  the 
number  of  frames  we  are  averaging.  Averaging  frames  serves  to  greatly  mitigate  the 
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ill  effects  of  speckle  noise  in  laser  illuminated  images  [42],  Finally,  averaging  frames 
also  serves  as  the  catalyst  for  using  the  average  model  for  the  atmosphere  when  trying 
to  enhance  the  spatial  and  temporal  resolution  of  the  3-D  image. 

An  additional  benefit  of  using  MCFA  images  is  that  it  will  also  serve  to  increase 
A4.  Given  the  variance  of  an  individual  image  as  shown  in  (2.13),  the  variance  of  a 
2-image  MCFA,  aNB 2,  will  be 

<JNB2  =  2K  ■  (2-15) 

If  we  can  now  assume  that  the  expected  photons  received  for  each  image  is  constant, 

Ki  =  K 2  =  K,  (2.16) 

we  can  perform  a  substitution  of  variables  to  show  that  A4  doubles  for  a  2-frame 
MCFA.  The  2-frame  MCFA  will  also  have  a  negative  binomial  distribution  with 

* f'(1  +  li)=iK1  +  ^)'  (217) 

where 

k'  =  K\  +  J?2  =  2k.  (2.18) 

Since  the  shape  of  the  negative  binomial  PMF  is  primarily  driven  by  A4,  this  increase 
in  A4  will  ultimately  make  our  Poisson  approximation  more  accurate. 

2.3  Experimental  Data  Collection  Systems 

This  section  will  present  the  two  camera  configurations  that  were  used  to  obtain 
experimental  results  to  support  this  dissertation.  All  of  the  equipment  was  com¬ 
mercially  available,  but  in  some  cases  modified  slightly  to  perform  in  accordance  the 
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expectations  designed  into  the  algorithms. 


2.3.1  2-D  Visible  Light  System. 

The  optical  configuration  shown  in  Figure  2.3  was  used  to  obtain  properly  sampled 
2-D  images.  The  experimental  results  obtained  with  this  setup  were  primarily  used  to 
verify  initial  concepts  on  the  impact  of  turbulence  on  images,  the  ability  to  measure 
r0  and  in  support  of  the  research  presented  in  Chapter  III.  The  specifications  for  this 
setup  are  listed  in  Table  2.1.  The  camera  used  in  this  configuration  allowed  for  16-bit 


Figure  2.3:  Experimental  sensor  setup  consists  of  a  Celestron®  NexStar®  6SE  1.5 
m  focal  length  telescope  with  a  mask  to  reduce  the  aperture  to  5  cm,  and  an  Orion® 
Starshoot™  G3  monochrome  camera. 


Analog  to  Digital  (A/D)  conversion  and  experiments  show  that  it  acts  as  a  photon 
counting  device  in  lower  intensity  regions.  This  allows  for  utilization  of  the  Poisson 
and  negative  binomial  model  assumptions  without  applying  a  conversion  factor  to  the 
digitized  images.  However,  as  the  detector  approached  higher  intensity  thresholds,  the 
conversion  between  digital  counts  and  photons  became  non-linear.  For  that  reason, 
images  were  taken  in  low  light  conditions  to  ensure  the  various  techniques  developed 
through  this  research  could  be  applied. 
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Table  2.1:  Optical  System  Specifications. 


Parameter  Name 

Defined  Value 

Mean  wavelength  (A) 
Detector  array  size 
Pixel  size 

Sensor  focal  length  (/) 
Aperture  diameter  ( D ) 

600  nm 

582  x  582 

8.3  /ini  x  8.6  /mi 
1.5  m 

5  cm 

2.3.2  3-D  FLASH  LADAR  System. 

The  3-D  FLASH  LADAR  sensor  used  for  this  research  was  a  modified  ASC 
Portable  3-D  FLASH  LADAR  Camera  Kit™.  The  sensor  used  all  of  the  standard 
components  in  this  commercially  available  sensor;  however,  they  were  oriented  into  a 
different  configuration  for  the  USAF  Test  Pilot  School  (TPS)  as  shown  in  Figure  2.4. 
The  system  specifications  as  configured  are  listed  in  Table  2.2. 


(a)  Standard  camera  configuration  (b)  TPS  camera  configuration 

Figure  2.4:  (a)  This  is  the  standard  configuration  of  the  camera  available  from  ASC. 

(b)  In  order  to  meet  USAF  TPS  requirements,  the  camera  components  were  oriented 
into  a  brassboard  design  to  allow  for  installation  into  an  airborne  pod. 

The  experimental  3-D  FLASH  LADAR  sensor  possessed  a  fixed-focus  lens,  with 
focus  set  to  infinity.  While  this  configuration  would  be  advantageous  for  installation 
into  an  airborne  pod  where  long  slant  ranges  to  the  target  are  common,  it  did  create 
some  challenges  in  trying  to  establish  suitable  experimental  target  configurations. 


25 


Table  2.2:  ASC  Portable  3-D  FLASH  LADAR  System™Specifications. 


Known  System  Parameters 

Parameter  Name 

Defined  Value 

Frames  per  image 

20 

Mean  laser  wavelength  (A) 

1.57  /jrn 

Measured  sample  rate 

434.5  MHz 

Energy  per  pulse  (Et) 

0.025  J 

Sensor  pulse  width 

4.7x10-9s 

Detector  array  size 

128  x  128 

Pixel  size 

100  /jrn  x  100  /jrn 

Lens  Parameters 

Parameter  Name 

Defined  Value 

Sensor  focal  length  ( fi ) 

250  mm 

Aperture  diameter  ( D ) 

12  cm 

Instantaneous  Field  of  View  (iFOV) 

3° 

Each  individual  3-D  image  consisted  of  20  2-D  frames.  In  SULAR  mode  where 
each  pixel  within  the  detector  is  set  to  a  fixed  range  gate,  this  would  result  in  a 
significant  limitation.  Based  on  the  sample  rate,  the  range  gate  would  be  limited  to 
approximately  6.56  m  with  each  2-D  frame  separated  by  0.345  m.  Fortunately  the 
camera  has  a  mode  of  operation  known  as  “STOP”  mode.  In  this  mode  of  operation, 
3-D  images,  each  with  20  frames,  can  be  taken  with  various  starting  points  for  the 
range  gate.  The  3-D  image  for  each  subsequent  STOP  will  be  formed  from  a  separate 
laser  pulse  and  will  provide  an  additional  8  frames  of  data  extending  the  range  gate 
by  approximately  2.76  m  per  STOP.  Figure  2.5  illustrates  how  the  composite  image 
is  reassembled  to  achieve  a  longer  range  gate.  The  camera  has  been  tested  with  a 
total  of  9  additional  STOPs  for  a  total  range  gate  of  31  m.  In  this  mode  of  operation 
we  are  truly  only  limited  in  the  sensor’s  ability  to  quickly  read  and  store  a  significant 
amount  of  digitized  data.  While  this  research  appears  somewhat  constrained  by  the 
requirement  of  using  the  SULAR  mode  of  operation,  techniques  such  as  STOP  mode 
mitigate  this  constraint. 
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Figure  2.5:  Assembly  of  Composite  3-D  Image  Using  Two  Additional  STOPS. 

2.3.3  Speckle  Parameter  Estimation  and  Photon  Calibration. 

The  algorithms  developed  in  this  research  are  based  on  statistical  models  for  the 
photon  arrival  rates  of  partially  coherent  and  incoherent  light.  Therefore,  a  calibration 
of  both  experimental  sensors  was  warranted.  For  instance,  the  iterative  algorithms 
will  use  a  stopping  criteria  that  depends  on  a  comparison  of  the  variance  between 
the  non-noisy  estimates  and  the  measured  data.  Additionally,  one  of  the  algorithms 
depends  on  the  ability  to  calculate  likelihood  across  a  range  of  parameters.  Therefore, 
we  must  understand  the  sources  of  variation  in  the  measured  data,  and  the  overall 
number  of  photons  received. 

Digital  cameras  commonly  produce  images  by  computing  a  digital  count  that 
corresponds  with  the  photons  received  at  the  detector  plane  over  the  integration 
period.  Unless  the  camera  is  a  true  photon  counting  device  where  the  digital  count 
is  equivalent  to  the  photons  received,  the  statistics  will  be  skewed.  The  variance  of  a 
random  variable  X  is  defined  as 

°x  =  E{{. X-Hx)2],  (2.19) 

where  fix  is  the  mean  of  the  random  variable  X.  Modeling  the  received  photons  with 
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the  random  variable  Y,  a  scaling  by  some  factor,  c,  will  yield 


Y  =  cX.  (2.20) 

In  converting  to  camera  counts,  the  variance  of  this  scaled  value,  cry,  will  become 

a2Y  =  E  [(Y  -  /iy)2]  =  E  [{cX  -  cnxf] 

4  =  (c2)  E  [(X  -  iiX)2]  (2.21) 

a'y  =  c2a2x. 

By  imposing  a  scaling  factor  when  converting  to  camera  counts,  the  variance  too  is 
scaled  by  the  square  of  that  scaling  factor.  Therefore,  if  the  distribution  of  the  light 
received  should  be  Poisson,  a  scaled  representation  of  the  photons  received  in  the 
form  of  camera  counts  will  clearly  not  be  Poisson. 

For  a  camera  that  is  measuring  natural  light,  the  conversion  of  camera  counts  is 
rather  simple.  Since  a  Poisson  process  will  have  a  mean  equal  to  the  variance,  we 
only  need  to  find  a  scaling  factor  that  allows  for  this  condition  to  hold.  However, 
caution  must  be  exercised  to  ensure  that  conversion  is  linear  for  the  intensity  range 
of  interest,  or  the  non-linearity  will  need  to  be  addressed.  For  the  experimental 
configuration  presented  in  Figure  2.3,  this  calibration  was  accomplished  by  first  taking 
a  series  of  images  with  the  lens  cap  on  in  order  to  capture  the  variance  associated 
with  the  detector  bias.  A  series  of  100  short  exposure  images  was  then  captured 
for  a  uniformly  illuminated  step  target.  This  step  target  captured  the  minimum 
and  maximum  intensity  range  that  we  expected  in  our  experiments.  The  mean  and 
variance  were  then  computed  for  this  series  of  images  in  an  attempt  to  find  the  scaling 
factor  as  shown  in  Figure  2.6.  Fortunately  the  mean  and  variance  were  found  to  be 
approximately  equal  for  the  intensity  ranges  we  were  concerned  with.  Given  that  the 
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Figure  2.6:  (a)  Shows  the  step  target  used  to  capture  the  minimum  and  maximum 

expected  intensity,  (b)  Comparison  of  pixel  mean  and  variance  for  a  1-D  slice  through 
the  middle  of  the  image. 


camera  has  a  16-bit  A/D  converter,  intensity  values  can  range  from  0  -  65535.  Our 
experimentation  found  that  as  long  as  an  individual  pixel  intensity  remained  below 
2,000  camera  counts,  there  was  approximately  a  1  to  1  ratio  for  camera  counts  to 
photons. 

With  the  calibration  for  natural  light  achieved,  the  calibration  for  the  coherence 
parameter,  Ad,  of  the  laser  illuminator  used  in  the  2-D  experiments  was  possible. 
Given  that  the  partially  coherent  light  will  have  a  negative  binomial  distribution 
we  simply  needed  to  find  the  Ad  that  minimized  the  Mean  Square  Error  (MSE) 
between  the  variance  computed  from  (2.13)  and  the  measured  variance.  For  the  2-D 
experiments,  the  remote  scene  was  illuminated  using  a  laser  with  a  wavelength  of 
630nm  and  a  measured  coherence  parameter,  Ad  =  10. 

The  calibration  of  the  3-D  FLASH  LADAR  sensor  was  more  complex.  Previous 
research  has  relied  on  a  two  step  calibration  process  that  was  similar  to  what  was  used 
for  the  2-D  experimental  setup  [47].  However,  this  technique  depends  on  the  ability 
to  collect  images  that  are  illuminated  by  incoherent  light.  This  was  accomplished  by 
pointing  the  sensor  at  a  bright  scene,  and  collecting  images  with  the  laser  illumina- 


29 


tor  inhibited.  Unfortunately,  this  technique  was  not  possible  with  the  3-D  FLASH 
LADAR  camera  used  for  this  research.  It  is  unknown  if  the  lens  on  this  particular 
setup  used  a  filter  to  attenuate  wavelengths  outside  of  a  narrow  band  around  the  laser 
illuminator,  as  attempts  to  calibrate  the  photon  to  camera  count  ratio  with  natural 
light  were  inconsistent.  Therefore,  a  scheme  was  required  that  allowed  for  simultane¬ 
ous  calibration  of  the  camera  count  to  photon  calibration  and  coherence  parameter. 
A  technique  was  developed  that  uses  the  y2  goodness  of  fit  test  to  allow  for  this  cal¬ 
ibration.  Finally,  the  3-D  FLASH  LADAR  camera  has  numerous  different  detector 
sensitivity  settings  that  will  allow  the  camera  count  to  photon  count  calibration  to 
vary.  Previous  work  with  an  ASC  LADAR  system  [7],  and  information  provided  by 
the  manufacturer  led  us  to  believe  that  this  system  was  a  photon  counting  device. 
In  other  words,  we  expected  the  photon  to  camera  count  conversion  factor  to  be  1. 
However,  we  were  looking  for  a  technique  to  confirm  this  information  to  be  true  for 
the  system  used  for  this  research. 

Given  that  we  can  easily  collect  dark  images  with  the  3-D  FLASH  LADAR  camera 
by  placing  a  cover  over  the  lens,  it  was  possible  to  characterize  the  variance  attributed 
to  the  combined  camera  bias.  The  remaining  variance  in  the  images  can  primarily  be 
attributed  to  the  collection  of  the  partially  coherent  illumination.  In  regions  of  high 
intensity  returns  within  the  image,  the  variance  will  be  dominated  by  speckle  noise 
with  a  negative  binomial  distribution.  Additionally,  it  is  known  that  the  shape  of 
the  distribution  is  primarily  determined  by  the  coherence  parameter  [21].  Therefore, 
a  technique  was  derived  using  the  chi-squared  goodness  of  fit,  or  y2  test,  under  the 
hypothesis  that  we  could  use  this  test  to  first  identify  a  likely  candidate  for  the 
coherence  parameter  given  the  shape  of  the  distribution.  The  y2  goodness  of  fit  test 
statistic,  X2,  is 


x2  =  X 


(Oz  -  Ezf 
Ez 


(2.22) 
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where  we  form  an  expected  distribution  with  Z  bins  and  Oz  is  the  observed  number 
of  occurrences  and  Ez  is  the  expected  number  of  occurrences  for  a  particular  bin. 

Knowing  that  the  reflected  illumination  is  partially  coherent,  we  seek  to  find  out  if 
the  received  data  are  from  a  particular  negative  binomial  distribution.  We  still  have 
two  unknowns,  the  photon  scaling  parameter  and  Af.  However,  since  the  shape  of 
the  distribution  is  primarily  determined  by  At,  we  expected  that  the  y2  test  would 
reveal  similar  results  across  a  range  of  scaling  factors.  We  could  then  use  the  results 
for  At  to  directly  calculate  the  scaling  based  on  a  comparison  with  the  theoretical 
variance  in  (2.13).  If  our  X2  statistic  is  high,  we  can  conclude  that  the  data  are  not  a 
part  of  the  tested  distribution  and  accept  the  alternate  hypothesis,  Hi,  which  states 
we  have  a  lack  of  fit.  However,  if  our  X2  statistic  is  low,  we  fail  to  reject  the  null 
hypothesis,  H0.  Therefore,  we  are  looking  for  the  non-rejection  region  where  we  fail 
to  reject  H0  in  order  to  identify  the  most  likely  value  of  Ah  In  using  the  y2  test  to 
find  the  At  that  produces  the  best  fit,  we  test  two  hypotheses 


X2  <  y2  (1  —  cy  DOF)  :  Conclude  H0 
X2  >  y2  (1  —  a;  DOF)  :  Conclude  Hi, 


(2.23) 


with  a  level  of  significance,  a,  and  the  degrees  of  freedom  determined  by  the  number 
of  bins  selected  for  the  distribution  [36]. 

We  first  use  (2.12)  to  find  the  estimated  PMF  for  a  given  value  of  K  and  At.  We 
can  then  form  an  experimental  PMF  by  first  scaling  and  then  taking  a  histogram  of 
the  measured  data.  Both  the  scaling  factor  and  At  will  be  varied  over  a  predeter¬ 
mined  range.  Our  research  has  shown  that  the  change  in  the  non-rejection  region 
with  increasing  scale  factor  values  is  similar  to  a  decaying  exponential.  For  scale 
factors  below  the  true  value  we  get  high  predictions  for  At.  However,  as  the  scale 
factor  approaches  and  then  exceeds  the  true  value,  the  non-rejection  region  begins 
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to  stabilize.  Upon  conducting  numerous  simulations,  our  hypothesis  appears  to  be 
confirmed  in  that  the  true  value  of  Ai  is  generally  centered  within  the  non-rejection 
region  once  it  stabilizes.  This  allows  us  to  narrow  down  the  potential  combinations 
of  scale  factor  and  coherence  parameter  significantly. 

It  should  be  noted  at  this  point,  that  the  primary  goal  of  this  calibration  proce¬ 
dure  is  to  get  reasonably  close  to  the  true  values  of  photon  scaling  factor  and  Ai. 
Additionally,  we  wanted  to  find  an  acceptable  combination  of  scale  factor  and  Ai  that 
made  the  collected  data  appear  to  have  a  negative  binomial  distribution.  Depending 
on  the  scale  factor,  some  deviation  from  the  truth  in  the  estimates  with  this  technique 
would  be  expected  due  to  the  error  associated  with  scaling  also  known  as  quantization 
error.  Further,  as  observed  in  Figure  2.2,  as  Ai  grows,  minor  deviations  in  Ai  have 
little  impact  on  the  shape  of  the  distribution.  Experimental  results  demonstrate  that 
minor  deviations  in  both  the  photon  scaling  factor  and  Ai  have  minimal  impact  on 
the  final  results  produced  by  the  algorithms  presented  in  this  dissertation.  Once  the 
likely  region  for  Ai  is  identified,  we  can  compare  the  scaled  variance  for  the  measured 
data  according  to  (2.21)  with  the  theoretical  variance  according  to  (2.13). 

The  following  is  a  demonstration  of  the  aforementioned  calibration  technique  on 
simulated  data.  In  each  of  the  three  simulations,  1000  independent  samples  with 
a  negative  binomial  distribution  were  generated.  The  data  were  then  divided  by 
the  scale  factor  to  represent  simulated  camera  count  data.  Table  2.3  presents  the 
true  and  estimated  values  for  scaling  and  Ai  for  each  of  the  three  trials.  Figure  2.7 


Table  2.3:  Scaling  and  Coherence  Calibration  Simulation. 


True  Parameters  vs.  Estimates 

Trial  # 

Scale  Factor 

Coherence  Parameter 

True 

Estimate 

True 

Estimate 

1 

3 

3 

150 

150 

2 

2 

2.1 

250 

235 

3 

5 

5.2 

200 

200 
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presents  the  findings  from  each  of  the  primary  steps  listed  above.  For  each  of  the 
trials,  the  hypothesis  test  was  conducted  for  a  range  of  scale  factors  from  1-20  and 
a  range  of  A4  values  from  1  -  500.  The  true  value  for  A4  was  then  chosen  from  the 
center  of  the  stabilized  non-rejection  region.  Using  this  value  for  A4  we  also  show 
the  point  of  intersection  between  the  theoretical  variance  and  the  measured  variance 
for  observations  over  a  range  of  scale  factors.  The  estimated  scale  factor  to  convert 
camera  counts  to  photons  is  the  point  of  intersection. 
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Figure  2.7:  In  this  demonstration,  the  regions  in  red  are  where  the  null  hypothesis  can 
be  rejected.  The  regions  in  blue  are  where  we  fail  to  reject  the  null  hypothesis.  The 
white  box  indicates  the  scale  factor  and  M.  from  which  the  results  are  based.  The  plots 
on  the  right  hand  side  show  the  intersection  of  the  measured  variance  and  theoretical 
variance  for  a  range  of  scale  factors. 
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Experimental  measurements  generally  consisted  of  144  3-D  images.  Based  on 
the  simulations  presented  in  Figure  2.7  and  the  limited  number  of  samples  available, 
we  expected  a  wider  non- rejection  region  as  shown  in  Figure  2.8(a).  However,  we 
were  able  to  find  a  reasonable  range  for  A4  that  had  negligible  impact  on  the  overall 
results.  Based  on  the  simulation  results,  we  expect  that  more  samples  might  enable 
a  tighter  bound  on  the  prediction.  Unfortunately,  the  time  required  to  collect  this 
data  would  likely  introduce  additional  inconsistencies  in  the  data  that  would  need 


to  be  addressed.  Based  on  the  results  in  Figure  2.8,  we  expect  that  the  camera  is 
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Figure  2.8:  Hypothesis  test  performed  on  data  collected  with  ASC  3-D  FLASH 
LADAR  sensor.  Data  exhibits  a  similar  trend  to  what  was  shown  in  simulation.  The 
regions  in  red  are  where  the  null  hypothesis  can  be  rejected,  (a)  The  non-rejection 
region  appears  to  stabilize  centered  on  a  coherence  parameter  of  195.  (b)  Given  a 

coherence  parameter  of  195,  the  scale  factor  is  approximately  1. 


a  photon  counting  device  with  a  A4  of  approximately  195.  Varying  the  value  of  M. 
from  100  to  300  will  vary  the  scale  factor  from  1.3  to  0.8  respectively.  However,  in  all 
cases  where  we  fail  to  reject  the  null  hypothesis,  the  collected  data  exhibit  negative 
binomial  characteristics  with  the  specified  scaling  and  coherence  parameter. 

As  a  final  note,  it  would  likely  be  preferential  to  calibrate  the  sensor  during  the 
design  process.  For  instance,  the  detector  could  be  calibrated  with  incoherent  il¬ 
lumination  and  a  separate  calibration  could  be  performed  on  the  coherence  of  the 
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laser  illuminator.  The  technique  described  above  using  the  y2  test  should  be  effec¬ 
tive  in  scenarios  where  prior  calibrations  are  not  available.  However,  caution  must 
be  exercised  to  avoid  distorting  the  variance  measurements  clue  to  imperfect  image 
registration.  Any  registration  errors  will  cause  a  spike  in  the  measured  variance  as 
shown  in  Figure  2.6  at  the  edge  of  the  step  target.  Therefore,  we  want  to  make  our 
calibration  measurements  near  largely  uniform  areas  within  the  target. 
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2.4  Iterative  Algorithm  Stopping  Criteria 


A  common  challenge  with  iterative  algorithms  and  a  topic  that  has  generated  a 
significant  amount  of  recent  research  involves  the  selection  of  a  stopping  criteria  [2], 
[12],  [19],  [56].  The  research  conducted  by  MacManus  and  presented  in  Chapter  III 
depended  primarily  on  the  stopping  criteria  to  arrive  at  an  estimate  for  Fried’s  Seeing 
Parameter,  tq.  In  the  multi-surface  ranging  algorithm  detailed  in  this  dissertation,  a 
key  feature  is  the  ability  to  maximize  likelihood  for  the  correct  value  of  r0.  However, 
since  likelihood  will  increase  with  each  iteration  even  if  the  wrong  value  of  ro  is 
selected,  we  may  find  situations  where  likelihood  is  higher  for  a  fixed  number  of 
iterations  with  the  wrong  r0  than  it  is  for  the  correct  value  of  r0.  This  likely  has  to 
do  with  the  rate  at  which  the  algorithm  converges.  Based  on  experience  with  the 
algorithms,  higher  values  of  r0  generally  allow  faster  convergence  than  lower  values. 
However,  if  the  algorithm  is  allowed  to  iterate  for  too  long,  we  can  distort  the  results. 

In  the  algorithms  employed  for  both  2-D  and  3-D  parameterized  blind  deconvo¬ 
lution,  the  model  for  the  received  image  considered  a  mean  intensity  convolved  with 
a  PSF  and  further  degraded  by  some  additive  amount  of  noise.  We  can  solve  for  the 
mean  level  of  this  noise.  However,  if  the  iterative  algorithms  are  not  stopped  at  the 
appropriate  time,  the  estimates  for  image  intensity  will  ultimately  be  distorted  by 
trying  to  fit  an  estimate  to  the  noise.  This  research  will  employ  a  Convergence  of 
Variance  (CoV)  technique  for  the  stopping  criteria.  This  technique  has  been  employed 
in  similar  research  in  the  past  with  success  [45],  [48].  The  data  model  is  provided  in 
Figure  2.1.  Given  this  relationship,  iterations  will  continue  until  the  MSE  between 
the  collected  data  and  the  image  estimate  is  lower  than  the  estimated  data  variance. 
For  a  3-D  image  this  criteria  can  be  represented  as 

K  M  K  M 

EE  (d(u,v,rk)  ~  I  ( u,v,  rk ))2  <  EE  V(u,v,rk),  (2.24) 

k=  1  u,v=l  k= 1  u,v=l 
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where  V  is  the  estimated  data  variance.  At  this  point  we  can  assume  that  any  further 
iterations  would  only  serve  to  fit  the  estimates  to  the  noise,  invalidating  the  model 
used. 

Multiple  techniques  exist  for  identifying  the  estimated  data  variance.  First  of  all, 
we  could  directly  compute  the  variance  if  we  have  a  series  of  images  taken  of  the 
same  scene.  Alternatively,  if  the  sensor  can  be  accurately  calibrated  and  we  have 
an  understanding  of  the  nature  of  the  light  being  detected,  we  can  mathematically 
predict  the  estimated  data  variance.  Given  a  partially  coherent  illumination  source, 
the  noise  in  the  images  will  likely  be  dominated  by  speckle  noise  where  the  image 
variance  can  be  accurately  modeled  with  the  negative  binomial  PDF  [58].  Therefore, 
the  expected  variance  for  each  image  pixel  and  range  slice  would  be 


VNB  (u,  v,  rk)  =  dNB  (u,  v,  rk) 


dNB(u,v,rk ) 

M 


(2.25) 


where 

dNB  (u,  v,  rk)  =  d  ( u ,  v,  rk)  -  B  (u,  v) .  (2.26) 

Further  refinement  can  be  given  to  this  model  when  the  variance  of  the  detector  bias 
is  also  considered.  In  FLASH  LADAR,  detector  bias  can  be  measured  by  taking 
images  without  firing  the  laser  and  is  generally  modeled  with  the  Poisson  PMF  [58]. 
Therefore  the  overall  variance  would  be 


V  (u,  v,  rk)  =  B  (u,  v )  +  VNB  ( u ,  v,  rk) .  (2.27) 

The  actual  distribution  of  the  variance  will  be  sensor  dependent;  however,  the  model 
in  (2.27)  accurately  approximated  the  variance  of  the  calibrated  3-D  FLASH  LADAR 
sensor  used  for  the  experimentation  in  this  research. 
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2.5  Seeing  Parameter  Measurement  from  Collected  Imagery 

In  Figure  2.9  we  demonstrate  the  ability  to  measure  r0  by  measuring  the  step 
response  from  the  collected  image.  The  impulse  response  is  then  found  by  taking  the 
derivative  of  this  measured  step  response.  Once  we  have  the  impulse  response  we 
can  vary  Tq  per  the  relationship  in  (2.9)  to  find  the  theoretical  total  OTF  that  min¬ 
imizes  the  error  between  the  measured  impulse  response  and  the  theoretical  impulse 
response.  This  demonstration  was  accomplished  by  blurring  a  perfect  step  target 
with  an  average  short  exposure  OTF  using  an  ro  of  0.0015m. 
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Figure  2.9:  (a)  Step  target  without  any  blur,  (b)  Step  target  blurred  using  a  short 

exposure  OTF  with  an  ro  of  0.0015m.  (c)  Measured  step  response  by  looking  at  the 
horizontal  change  in  intensity,  (d)  Impulse  response  computed  by  taking  the  first 
derivative  of  the  step  response,  (e)  Frequency  response  computed  by  taking  the  Fourier 
transform  of  the  impulse  response. 
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2.6  Previous  Research 


The  following  section  will  summarize  and  contrast  previously  published  research 
on  topics  pertaining  to  this  dissertation.  On  the  specific  topic  of  parameterized  blind 
deconvolution  with  multi-surface  ranging,  no  other  research  could  be  found.  However, 
this  dissertation  levied  ideas  from  several  sub-topics  in  the  final  formulation  of  the 
overarching  research  goal. 

2.6.1  Blind  Deconvolution. 

The  importance  of  image  deblurring  is  evident  based  on  the  wealth  of  research 
conducted  on  the  topic  of  blind  deconvolution.  Due  to  the  ill-posed  nature  of  blind 
deconvolution  with  2-D  images,  it  is  mathematically  impossible  to  directly  solve  for 
the  PSF  impacting  collected  images  when  noise  is  present  [38].  Despite  this  hurdle, 
numerous  algorithms  have  been  developed  to  circumvent  these  mathematical  chal¬ 
lenges  with  considerable  success  by  making  various  assumptions  or  approximations 
[5],  [11],  [25],  [35],  [37],  [61],  [60],  [64]  and  [68].  While  the  problem  can  be  extremely 
challenging  with  2-D  images,  Millane  et  al.  realized  that  working  with  3-D  images 
presented  the  potential  to  reduce  common  constraints  such  as  sensor  sampling  re¬ 
quirements  and  convergence  time  on  iterative  algorithms  [49]. 

Complexity  of  the  algorithm,  applicability  to  certain  classes  of  blurring  functions 
and  the  time  required  to  perform  the  necessary  operations  are  common  concerns 
associated  with  blind  deconvolution  algorithms.  Kundur  and  Hatzinakos  authored 
a  pair  of  articles  that  provide  an  exceptional  summary  of  the  challenges  associated 
with  blind  deconvolution  and  the  various  strategies  commonly  employed  [33],  [34], 
The  articles  highlighted  numerous  considerations  for  image  processing  applications 
that  will  directly  affect  this  research.  First,  due  to  the  ill-conditioned  nature  of  the 
problem,  small  perturbations  in  the  received  data  can  lead  to  large  deviations  in  the 
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estimates  produced  by  various  algorithms.  As  an  extension  to  this  concept,  noise 
amplification  is  especially  likely  [15].  This  research  reduces  the  issues  associated 
with  noise  amplification  through  selection  of  an  iterative  algorithm  stopping  criterion 
based  on  noise  variance,  and  application  of  a  constraint  on  the  amplitude.  Another 
common  issue  with  blind  deconvolution  is  that  multiple  solutions  are  likely  to  exist 
[4],  [29].  This  research  will  minimize  the  chances  for  convergence  to  an  undesirable 
solution  through  educated  initialization  and  applied  constraints. 

In  many  remote  sensing  applications,  processing  time  and  complexity  are  the 
primary  design  considerations  for  an  algorithm.  As  an  example,  the  APEX  method 
is  a  non-iterative,  blind  deconvolution  technique  that  can  enhance  certain  classes  of 
imagery  in  near  real  time  [8],  [9].  The  technique  operates  on  a  restricted  class  of 
blurs,  in  the  form  of  2-D  radially  symmetric,  bell-shaped,  heavy-tailed  probability 
density  functions.  While  this  technique  is  very  effective  on  certain  classes  of  imagery, 
it  does  not  provide  an  estimate  of  r0  as  given  in  the  CoV  technique.  A  comparison  of 
this  work  to  the  CoV  technique  was  provided  in  [45].  As  previously  mentioned,  the 
ability  to  recover  ro  has  applicability  beyond  the  primary  focus  of  conducting  blind 
deconvolution. 

The  work  reported  on  in  [10]  and  [17]  is  significant  because  it  shows  the  value  of 
the  Richardson- Lucy  (RL)  Filter  for  iteratively  deblurring  images.  Variations  on  this 
technique  are  employed  throughout  this  research  in  conjunction  with  OTFs  param¬ 
eterized  by  ro-  One  technique  for  simplifying  this  problem  considers  a  parameteri¬ 
zation  of  the  OTF  [41],  [43]  and  [45].  This  work  will  extend  the  concepts  presented 
by  MacDonald  and  MacManus  to  show  that  working  with  3-D  images  presents  new 
opportunities  that  were  previously  mathematically  ill-posed  as  theorized  by  Millane 
[49], 

The  parameterized  OTF  is  chosen  to  simplify  the  structure  of  the  unknown  blur- 
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ring  function  to  a  certain  class  of  atmospheric  models.  The  goal  of  this  simplification 
is  to  make  the  algorithm  more  suitable  for  implementation  in  a  tactical  environment 
where  near-real-time  operation  is  required.  When  using  MCFA  compilations  of  im¬ 
ages  where  each  individual  image  has  an  exposure  time  of  less  than  jA  of  a  second, 
the  average  short-exposure  OTF,  Hse,  is  reduced  to  a  function  of  a  single  unknown 
parameter,  Fried’s  seeing  parameter,  r0,  as  shown  in  (2.10).  A  similar  average  OTF 
has  been  discussed  for  the  long  exposure  case  (2.11)  as  well.  While  there  are  benefits 
to  using  long-exposure  imaging  in  certain  scenarios  such  as  astrophotography,  the  loss 
of  frequency  content  is  often  an  undesired  side-effect.  In  a  tactical  military  scenario 
or  any  other  dynamic  environment,  it  would  be  impractical  if  not  impossible  to  point 
a  sensor  at  a  target  long  enough  to  warrant  the  use  of  long-exposure  imaging. 

Figure  2.10  compares  the  frequency  response  that  is  diffraction  limited,  to  the 
frequency  response  considering  a  long  and  short-exposure  OTF  for  two  levels  of  at¬ 
mospheric  seeing.  It  is  evident  from  a  frequency  content  standpoint  that  there  are 
inherent  benefits  with  using  properly  registered  short-exposure  images.  The  param¬ 
eters  for  this  demonstration  were  chosen  to  match  those  in  Table  2.1  that  will  be 
employed  to  obtain  the  experimental  results  shown  in  Chapter  III.  As  expected,  the 
short-exposure  OTF  is  very  close  to  the  diffraction  limited  OTF  when  ro  =  5  cm. 
However,  the  long-exposure  OTF  reveals  a  significant  attenuation  in  high  frequency 
content.  Higher  levels  of  turbulence  yield  a  higher  loss  in  frequency  content  for  the 
long-exposure  scenario,  and  significant  attenuation  of  high  frequency  content  for  the 
short-exposure  scenario  as  shown  in  Figure  2.10(b). 

One  of  the  primary  benefits  associated  with  parameterized  blind  deconvolution  is 
computational  efficiency.  While  the  most  significant  improvement  is  realized  through 
deconvolution  with  the  actual  instantaneous  OTF,  simultaneously  computing  the  in¬ 
stantaneous  OTF  and  deblurring  an  image  is  computationally  intensive.  However, 
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(a)  OTF  Comparison  for  r0=5  cm  (b)  OTF  Comparison  for  r0= 2  cm 


Figure  2.10:  (a)  Comparison  of  the  frequency  response  for  this  sensor  given  a  D/r$ 

ratio  of  1.  (b)  Comparison  of  the  frequency  response  for  this  sensor  given  a  D/ro  ratio 
of  2.5. 

this  research  will  show  that  deconvolution  with  an  average  theoretical  model  for  the 
OTF  can  produce  significant  gains  in  image  quality  while  reducing  the  computational 
burden.  Additionally,  in  scenarios  where  the  average  atmosphere  remains  relatively 
constant,  the  problem  can  almost  be  treated  as  a  deconvolution  problem  rather  than  a 
blind  deconvolution  once  the  initial  parameterized  blind  deconvolution  is  completed. 


2.6.2  Multiple  Surface  Ranging. 

Variations  of  LADAR  technology  have  proven  useful  in  a  myriad  of  civilian  and 
military  applications.  Common  civilian  uses  are  3-D  mapping  of  surfaces,  autonomous 
vehicle  navigation  and  forestry  classification.  Military  applications  include  tasks  such 
as  targeting  and  autonomous  aerial  refueling.  Current  military  application  may  be 
limited  clue  to  sensor  capability;  however,  as  the  technology  improves,  additional 
emphasis  will  be  given  to  incorporating  LADAR  onto  new  or  existing  platforms. 
Previously  conducted  research  and  the  demonstrated  ranging  accuracy  of  LADAR 
systems  have  forged  an  interest  in  a  myriad  of  defense  applications. 

Currently  there  is  a  significant  desire  within  the  Department  of  Defense  (DoD)  and 
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DHS  to  employ  a  3-D  imaging  sensor.  Provided  spatial  resolution  can  be  improved  to 
be  comparable  with  currently  employed  passive  sensors,  3-D  FLASH  LADAR  technol¬ 
ogy  is  a  likely  successor.  One  of  the  primary  motivators  for  this  move  is  that  previous 
simulations  and  experiments  with  both  ALS  and  FLASH  systems  have  demonstrated 
the  ability  to  image  through  foliage  canopies,  camouflage  netting  or  various  other 
obscurations. 

In  LADAR  imaging,  a  common  method  for  obtaining  the  ranges  to  multiple  sur¬ 
faces  per  detector  is  to  fit  a  Gaussian  mixture  (2.2)  to  the  received  pulse  through  the 
use  of  various  techniques  such  as  an  EM  algorithm  [6],  [27],  [55]  and  [67].  Hernandez 
et  al.  developed  a  multiple  surface  ranging  algorithm  that  relied  on  reversible  jump 
Markov  chain  Monte  Carlo  techniques  [24],  While  their  work  shows  great  promise 
for  many  applications,  the  complexity  of  the  algorithm  does  not  lend  itself  to  near 
real  time  image  processing.  Additionally,  most  algorithms  operate  on  each  pixel  in 
isolation  and  do  not  consider  the  spatial  or  temporal  effects  from  interaction  with 
neighboring  surfaces.  While  the  techniques  may  give  additional  insight  to  candidate 
ranges  for  an  image,  depending  on  the  severity  of  the  diffraction  effects,  the  surfaces 
visible  in  the  collected  data  may  vary  significantly  from  reality.  The  primary  goal  of 
this  portion  of  the  research  was  to  develop  an  EM  solution  which  could  accurately 
discriminate  the  range  to  multiple  true  surfaces  per  pixel  while  discarding  false  re¬ 
turns  due  to  diffraction.  Since  we  are  considering  statistical  independence  for  each 
detector  in  the  APD  array,  and  the  received  pulse  is  a  composition  of  temporally 
displaced  Gaussians,  ideas  and  techniques  for  pulse  estimation  can  be  leveraged  from 
numerous  other  disciplines.  Gaussian  decomposition  has  applications  in  nearly  any 
case  where  the  collected  data  can  be  decomposed  into  numerous  subsets,  each  with 
an  approximately  normal  distribution. 

Fortunately  Gaussian  decomposition  or  mixture  modeling  is  studied  a  great  deal 
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in  the  literature,  but  numerous  challenges  are  commonly  cited  with  reference  to  de¬ 
veloping  an  estimator  for  the  Gaussian  mixture.  Expectation  Maximization  (EM)  is 
a  common  technique  employed  to  isolate  the  parameters  of  interest  in  the  received 
data  [24],  [46].  Zhuang  points  out  that  common  challenges  are  employing  estima¬ 
tors  when  the  number  of  components  in  the  mixture  are  unknown,  or  when  mixture 
components  may  merge  clue  to  their  individual  parameters  [69].  Vlassis  and  Nikas 
approached  the  first  problem  with  a  greedy  EM  technique.  They  performed  an  it¬ 
erative  process  where  the  number  of  components  was  incrementally  increased  until 
the  number  corresponding  with  the  solution  with  the  highest  likelihood  was  obtained 
[66].  With  3-D  FLASH  LADAR,  this  process  could  be  extremely  time  consuming  as 
the  number  of  components  grows,  therefore  making  this  an  impractical  solution  for 
this  work.  Rather,  this  research  will  propose  the  establishment  of  an  upper  bound, 
and  then  refining  the  estimate  based  on  pulse  parameter  estimates. 

The  approach  employed  for  multi-surface  ranging  is  similar  to  the  EM  approach 
derived  by  Dolce  for  fusing  2-D  and  3-D  LADAR  data  [16].  However,  the  work  is  not 
equivalent  because  Dolce’s  work  does  not  account  for  the  possibility  of  multiple  sur¬ 
faces  per  pixel,  nor  does  it  provide  an  estimate  for  the  pulse  amplitude  or  pulse  width. 
Dolce’s  algorithm  was  only  concerned  with  a  single  pulse  detector.  For  this  reason, 
he  was  able  to  use  the  amplitude  received  from  the  Richardson-Lucy  deconvolution 
of  the  2-D  image.  However,  since  we  would  like  to  solve  the  problem  where  multiple 
surface  returns  may  be  received  by  each  detector,  we  will  not  have  that  luxury.  An 
interesting  topic  for  future  research  would  be  to  combine  the  work  presented  in  this 
dissertation  with  the  work  conducted  by  Dolce. 

A  considerable  amount  of  previous  research  focused  on  LADAR  technology  has 
been  devoted  to  enhancing  the  capability  of  ALS  systems  such  as  in  [28],  [32],  [46] 
and  [63].  One  advantage  for  ALS  sensors,  is  that  their  spatial  resolution  is  primarily 
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determined  by  the  system’s  ability  to  finely  scan  over  the  target  area  and  accurately 
reconstruct  the  received  pulses  into  an  image.  Mallet  and  Bretar  provide  a  good  sum¬ 
mary  of  various  techniques  and  successes  that  have  been  realized  with  this  technology 
with  respect  to  multiple  surface  ranging  [46].  Unfortunately,  ALS  systems  commonly 
require  a  considerable  amount  of  time  to  form  an  image,  and  the  targeted  scene  must 
remain  fairly  constant  during  the  imaging  period.  In  addition  to  the  constraints  on 
the  targeting  area,  ALS  sensors  also  suffer  from  added  complexity  due  to  the  require¬ 
ment  for  a  tracking  and  stabilization  system.  As  pointed  out  by  Halmos,  one  of  the 
primary  driving  factors  for  the  interest  in  FLASH  systems  is  that  “size  can  be  dra¬ 
matically  reduced  by  eliminating  the  stabilization  subsystem  that  can  be  a  large  part 
of  the  LADAR  implementation  cost  and  size”  [23]. 

As  3-D  FLASH  technology  improves,  the  benefits  of  ALS  systems  without  the 
detractors  may  be  possible.  In  addition  to  the  research  tailored  specifically  to  ALS 
systems,  the  following  work  has  been  devoted  to  multiple  surface  ranging  with  FLASH 
or  hybrid  systems.  By  gathering  polarization  information  from  various  surfaces  in  con¬ 
junction  with  3-D  FLASH  LADAR  data,  Murray  demonstrated  that  multiple  surfaces 
could  be  discerned.  Additionally,  the  undesired  obscuration  could  be  discarded  based 
on  a  priori  knowledge  of  the  polarization  of  the  desired  surface  [50].  In  addition  to 
the  more  common  imaging  through  camouflage  applications,  Gelbart  et  al.  used  a 
FLASH  system  to  image  through  ocean  environments  to  demonstrate  the  ability  to 
detect  sea  mines  or  obstructions  that  may  impact  a  beach  landing.  Gelbart  ascer¬ 
tained  the  number  of  surfaces  in  the  received  data  by  counting  the  zero-crossings  of 
the  first  derivative  of  the  received  pulse  [18].  An  observation  in  common  with  both 
efforts  was  that  the  spatial  resolution  for  collected  images  was  poor  due  to  hardware 
limitations  and  complexity  of  the  detector  array.  As  a  proof-of-concept,  the  JIG¬ 
SAW  LADAR  sensor  developed  for  the  Defense  Advanced  Research  Projects  Agency 
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(DARPA)  used  a  hybrid  of  both  FLASH  and  ALS  technology.  By  scanning  with 
an  8x128  detector  array,  the  sensor  could  rapidly  obtain  a  304x256  3-D  image  [40]. 
While  the  spatial  resolution  of  this  system  was  an  improvement,  the  added  complex¬ 
ity  and  increase  in  hardware  associated  with  tracking  and  image  stabilization  would 
likely  limit  the  technique’s  utility  for  many  applications  where  size  is  an  important 
design  consideration. 

The  previously  mentioned  research  highlights  the  potential  advantages  of  employ¬ 
ing  3-D  LADAR  or  more  specifically  3-D  FLASH  LADAR  on  future  sensor  platforms. 
Common  among  each  of  the  previously  mentioned  efforts  is  the  lack  of  addressing  the 
effects  of  imaging  through  a  turbid  medium.  The  techniques  and  associated  systems 
have  demonstrated  the  ability  to  identify  multiple  surfaces  per  detector.  However,  no 
previously  developed  algorithm  has  been  identified  that  accounts  for  the  spatial  mix¬ 
ing  associated  with  the  lowering  of  the  spatial  frequency  cutoff,  and  how  it  impacts 
multi-surface  ranging.  As  sensors  improve  in  maximum  attainable  spatial  resolution, 
this  impact  will  only  become  more  pronounced. 

2.6.3  3-D  FLASH  LADAR  Image  Enhancement. 

Recent  work  conducted  at  Utah  State  LIniversity  focused  on  improving  the  ability 
to  perform  Automated  Target  Recognition  (ATR)  with  3-D  LADAR  images  [54], 
The  focus  of  this  research  was  to  use  post  processing  techniques  to  enhance  the 
resolution  of  surface  edges  through  multi-surface  ranging  techniques.  While  the  focus 
of  the  research  has  similar  goals,  there  were  several  underlying  assumptions  that 
reduce  the  applicability  to  the  specific  problems  being  addressed  in  this  dissertation. 
Perhaps  the  most  difficult  assumptions  to  overcome  are,  that  the  range,  intensity, 
and  pointing  direction  of  each  return  was  known  and  that  the  surfaces  were  assumed 
to  have  constant  reflectivity  and  a  Lambertian  bi-directional  reflectance  distribution 
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function. 


Likely  the  most  similar  research  to  what  is  presented  in  this  dissertation  was 
produced  by  McMahon  [47],  [48].  McMahon’s  work  was  the  only  previous  research 
that  could  be  found  with  a  focus  on  improving  or  restoring  the  images  collected 
through  blind  deconvolution.  He  found  that  by  using  the  range  diversity  present 
in  3-D  FLASH  LADAR  images,  he  could  simultaneously  solve  for  the  instantaneous 
OTF  and  an  enhanced  model  for  the  3-D  image.  Several  differences  exist  between  his 
research  and  what  is  presented  in  this  dissertation.  First  and  foremost,  McMahon’s 
work  only  accounts  for  a  single  surface  per  detector.  He  ultimately  proposes  the 
topic  of  multi-surface  ranging  as  future  research.  Additionally,  he  produces  estimates 
for  the  instantaneous  OTF  rather  than  the  parameterized  average  OTF.  While  this 
technique  has  application  where  post-processing  can  be  conducted,  it  is  likely  too 
complex  to  be  conducted  in  a  near  real  time  fashion. 

Based  on  the  similar  structure  of  the  iterative  algorithm  developed  by  McMa¬ 
hon  [48]  and  the  Multi-Surface  Including  Diffraction  (MSID)  algorithm  presented  in 
Chapter  IV,  we  can  compare  rough  estimates  on  the  computational  power  required 
for  each  iteration.  Likely  the  most  complex  mathematical  component  in  each  algo¬ 
rithm  is  the  computation  of  the  2-D  Fast  Fourier  Transform  (FFT)  or  2-D  inverse 
Fast  Fourier  Transform  (iFFT).  The  remaining  operations  are  simply  pixel-by-pixcl 
additions,  multiplications  and  divisions.  The  2-D  FFT  can  be  used  to  efficiently  per¬ 
form  convolution  by  multiplying  the  2-D  FFT  of  both  components  and  then  taking 
the  2-D  iFFT  of  the  result.  Similarly,  the  use  of  2-D  FFTs  and  iFFTs  can  also  be  used 
to  speed  up  the  correlation  operations  throughout  the  algorithms  by  also  taking  the 
complex  conjugate  of  one  of  the  two  components.  Given  the  solutions  for  the  pulse 
shape,  PSF,  gain  and  bias  derived  by  McMahon  [48],  each  iteration  will  take  three  2- 
D  iFFTs,  and  six  2-D  FFTs.  For  each  level  of  seeing  tested,  an  iteration  of  the  MSID 
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algorithm  can  be  completed  with  just  two  2-D  iFFTs,  and  two  2-D  FFTs.  This  results 
in  a  roughly  56%  decrease  in  processing  time  per  iteration  for  the  MSID  algorithm. 
As  we  will  discuss  in  Chapter  V,  parameterized  blind  deconvolution  in  conjunction 
with  the  MSID  algorithm  is  performed  in  a  search  routine.  However,  with  the  use  of 
parallel  processing  and  enough  available  processing  threads,  the  search  routine  could 
be  accomplished  in  roughly  the  same  time  that  a  single  level  of  atmospheric  seeing 
could  be  tested. 

A  final  contrast  with  the  research  conducted  by  McMahon  is  that  his  work  does 
not  produce  direct  estimates  for  all  of  the  pulse  parameters.  Rather,  he  produces 
a  refined  estimate  for  the  pulse,  and  then  employs  a  ranging  algorithm  separately. 
Despite  the  differences,  McMahon’s  work  established  a  solid  foundation  from  which 
the  research  in  this  dissertation  is  derived. 
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III.  Convergence  of  Variance  for  Seeing  Parameter 

Estimation 


In  this  chapter,  two  previously  developed  image  reconstruction  algorithms  are 
presented  that  will  remove  the  effect  of  atmospheric  turbulence  on  2-D  images.  While 
neither  of  the  algorithms  presented  in  this  chapter  are  novel,  this  chapter  will  consist  of 
an  expansion  upon  previous  findings  and  a  discussion  /  demonstration  of  the  potential 
challenges  with  applying  the  techniques  to  3-D  imagery.  The  primary  focus  of  this 
portion  of  research  was  to  identify  a  blind  deconvolution  technique  that  could  be 
employed  in  a  tactical  military  environment  where  both  time  and  computational 
power  are  limited.  Additionally,  the  following  techniques  have  application  in  the 
measurement  of  atmospheric  seeing  conditions.  In  a  blind  deconvolution  fashion,  the 
algorithms  simultaneously  compute  a  high  resolution  image  and  an  average  model  for 
the  atmospheric  blur  parameterized  by  Fried’s  seeing  parameter. 

Convergence  of  Variance  (CoV)  presented  in  Section  2.4  as  a  stopping  criteria  has 
been  incorporated  with  success  in  various  blind  deconvolution  algorithms  [44],  [48]. 
MacManus  later  demonstrated  that  CoV  by  itself  allowed  for  identification  of  the  best 
model  for  the  PSF  parameterized  by  r o  [45].  The  novelty  in  his  approach  was  that 
it  did  not  assume  a  prior  distribution  for  the  seeing  parameter,  rather  it  assessed 
the  convergence  of  the  image’s  variance  as  the  stopping  criteria  and  identification  of 
the  proper  seeing  parameter  from  a  range  of  candidate  values.  Experimental  results 
show  that  the  CoV  technique  allows  for  estimation  of  the  seeing  parameter  accurate 
to  within  0.5  cm  and  often  even  better  depending  on  the  signal  to  noise  ratio. 

This  chapter  is  organized  as  follows:  Section  3.1  presents  background  material  on 
image  deconvolution  and  the  challenges  associated  with  blind  deconvolution,  in  Sec¬ 
tion  3.2  the  two  techniques  will  be  compared  on  a  variety  of  artificially  blurred  images, 
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Section  3.3  demonstrates  the  value  of  the  CoV  technique  on  actual  collected  imagery 
and  Section  3.4  presents  considerations  and  challenges  associated  with  extending  this 
technique  to  3-D  imagery. 

3.1  Image  Deconvolution 

The  process  of  deconvolution  is  commonly  performed  on  collected  images,  i,  that 
are  degraded  by  a  PSF  and  noise.  The  goal  of  the  process  is  to  recover  the  true 
representation  of  the  remote  scene.  The  notation  for  the  image  model  in  this  chapter 
will  differ  slightly  from  Chapter  II.  In  this  chapter  we  are  primarily  concerned  with 
2-D  images.  In  the  2-D  case,  the  model  for  the  received  data,  d,  is  the  convolution  of 
the  remote  object,  o,  with  the  PSF,  h,  with  an  additive  amount  of  noise,  n,  as  shown 
in 

(o  (8)  h)  +  n  —  d.  (3.1) 

For  the  2-D  case,  we  are  simply  accounting  for  a  received  intensity  during  the  detec¬ 
tors  integration  time  rather  than  trying  to  account  for  the  entire  pulse  shape  reflected 
off  of  the  remote  scene.  As  previously  mentioned,  there  are  a  plethora  of  algorithms 
designed  to  aid  in  this  problem.  Perhaps  the  most  widely  accepted  or  recognized  algo¬ 
rithm  for  image  deconvolution  when  the  collected  images  follow  a  Poisson  distribution 
is  the  RL  algorithm. 

3.1.1  Richardson-Lucy  Deconvolution  Algorithm. 

One  of  the  key  benefits  of  the  RL  algorithm, 

5  Oc  y)new  =  °  (x >  y)old  (u-X,V~y\  r0),  (3.2) 
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where 


I  (u,v\r0)  =  E[d  (u,v)\r0\  =  6  (x,  y)  h  (u  -  x,  v  -  y  |r0)jj  (3.3) 

x  y 

for  a  wide  variety  of  imaging  applications  in  the  presence  of  Poisson  statistics,  is  that 
if  the  algorithm  converges,  it  will  converge  to  the  MLE  [39],  [62],  In  (3.2),  6  is  an 
estimate  for  o  and  /  is  the  expected  value  of  the  intensity  received  at  each  detector 
pixel  given  a  specific  PSF.  Since  (3.2)  is  an  iterative  algorithm  dependent  upon  previ¬ 
ous  estimates  for  the  image,  6  (x,  y)Qid,  we  must  provide  some  sort  of  initialization  to 
the  algorithm.  Typically,  the  initial  value  for  the  image  estimate  is  just  set  to  equal 
the  collected  data. 

A  commonly  cited  drawback  to  the  RL  algorithm  is  the  noise  amplification  that 
occurs  as  the  number  of  iterations  increases  [10],  [15].  Noise  amplification  is  a  com¬ 
mon  problem  with  iterative  ML  algorithms  where  the  algorithm  is  attempting  to  fit 
the  estimate  to  a  particular  distribution  as  closely  as  possible.  Therefore,  we  must 
stop  iterations  before  noise  amplification  occurs.  Clearly  this  could  be  accomplished 
interactively  by  the  user;  however,  for  an  automated  or  blind  routine,  any  required 
user  interaction  would  be  undesirable.  The  method  proposed  in  this  work  will  rely 
on  the  convergence  of  the  estimate  of  the  noise  power  and  the  predicted  variance  of 
the  collected  data  to  cease  iterations. 

3.1.2  Blind  Estimate  of  Seeing  via  Maximum  a  Priori  Technique. 

The  idea  of  using  the  RL  algorithm  in  a  blind  fashion  to  de-blur  an  image  was 
previously  presented  by  Fish  et  al.  [17].  However,  their  work  did  not  employ  the 
theoretical  models  for  the  long  and  short-exposure  transfer  functions  parameterized  by 
r0  [21].  Perhaps  the  most  similar  work  to  this  portion  of  the  research  was  accomplished 
by  MacDonald.  He  developed  a  blind  technique  that  was  iterative  in  nature  like  the 
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RL  algorithm,  yet  he  considered  a  priori  information  for  images  distorted  by  speckle 
noise  following  more  of  a  negative  binomial  distribution  [44],  MacDonald  considered 
a  priori  information  for  the  distribution  of  ro  in  hopes  of  maximizing  the  likelihood 
at  the  appropriate  level  of  seeing. 

Assuming  the  collected  image  can  be  approximated  with  the  Poisson  PMF,  and 
if  we  assume  independence  of  the  measurements  for  every  pixel  in  the  detector  array, 
we  can  state  the  joint  probability  of  the  observed  noisy  image,  i,  as 


p[I  =  d  (u,  v) ;  V  (u,  u)]  =  Y\  J][ 


I  (u,  V  |r0  g  I(u,v\r0  ) 

d  ( u ,  v)! 


(3.4) 


Ideally,  we  would  like  to  maximize  the  likelihood  or  log  likelihood 


L{r0)  =  [d(u,v)  In  [I  (u,v  |r0)]  -  I  (u,v  |r0)  -  In  [d  (u,  u)!]],  (3.5) 

U  V 

over  a  range  for  r0  to  identify  the  appropriate  PSF.  Unfortunately,  likelihood  contin¬ 
ually  increases  with  ro  as  illustrated  by  the  following  example  in  Figure  3.1.  While 
the  RL  algorithm  maximizes  likelihood  for  a  given  PSF,  the  ML  solution  over  a  range 
of  r0  values  does  not  necessarily  occur  when  the  correct  PSF  parameterized  by  r0  is 
chosen.  This  is  the  direct  problem  that  MacDonald  sought  to  solve  by  introducing  a 
priori  information  for  the  distribution  of  r0. 
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(a)  Original  Image 


(b)  Blurred  Image  (r0= 2.5  cm) 


(c)  Log  Likelihood 


Figure  3.1:  (a)  The  original  image  without  the  added  effects  of  atmospheric  blurring, 

(b)  The  image  that  would  be  received  by  a  sensor  with  the  specifications  listed  in  Table 
3.1  given  an  ro  of  2.5  cm.  (c)  Likelihood  as  a  function  of  ro.  Zooming  in  on  the  plot  of 
likelihood  vs.  rg  will  show  a  gradual  increase  in  likelihood  with  increasing  rg. 
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MacDonald  hypothesized  that  a  distribution  could  be  applied  to  the  value  for  r0 
based  on  the  intuitive  observation  that  the  seeing  is  seldom  extremely  better  than 
the  average  and  can  often  be  worse.  The  form  of  the  probability  density  function  for 
the  random  parameter  r0  was  assumed  to  be 


/.Ro  (ro) 


e-M2(r0/ravg) 

r  avg/M2 


(3.6) 


where  ravg  is  the  average  atmospheric  seeing,  and  M 2  is  the  number  of  pixels  in 
the  detector  array.  In  situations  where  ravg  is  unknown,  we  initialize  the  value  to  be 
equal  to  the  aperture  diameter  [41] .  We  can  then  execute  the  likelihood  maximization 
strategy  produced  by  MacDonald.  If  the  value  of  r0  that  maximizes  likelihood  is  less 
than  the  initialized  value  for  ravg ,  we  set  the  value  of  ravg  to  the  new  r0  estimate 
and  recompute  the  estimate  for  ro-  We  continue  this  process  until  the  estimate  for 
r0  is  equal  to  ravg.  When  applying  this  technique  to  the  example  illustrated  in 
Figure  3.1,  we  see  that  likelihood  is  maximized  near  the  correct  value  of  r o  as  shown 
in  Figure  3.2.  While  the  technique  is  successful  in  this  scenario,  two  mathematical 


(a)  Prior  for  r0  (b)  Log  Likelihood  Using  the  Prior  for  r0 


Figure  3.2:  (a)  Log  likelihood  of  the  exponential  prior  as  a  function  of  ro.  (b)  Overall 

log  likelihood  with  the  addition  of  the  prior.  With  the  addition  of  the  prior,  likelihood 
is  maximized  for  the  correct  value  of  rp. 
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challenges  remain.  First,  the  choice  of  an  exponential  distribution  for  r0  is  probably 
inaccurate.  Extremely  low  values  of  ro  are  expected  to  be  nearly  as  unlikely  as  high 
values.  Second,  the  effect  of  scaling  the  exponential  density  function  by  the  number 
of  pixels  in  the  detector  array  is  difficult  to  justify  for  partially  illuminated  scenes 
such  as  astronomical  images.  The  CoV  technique  will  remove  the  requirement  for  the 
prior  distribution  on  r0  and  allow  us  to  directly  converge  to  the  correct  value. 

3.1.3  Blind  Estimate  of  Seeing  via  CoV  Technique. 

The  CoV  technique  works  on  the  premise  of  searching  for  the  best  possible  PSF 
parameterized  by  r0  in  the  amount  of  time  available  for  processing.  Given  more  time, 
the  technique  will  provide  a  more  refined  estimate  for  rp.  We  will  first  explain  this 
technique  in  more  detail  using  the  assumption  that  the  images  collected  follow  the 
Poisson  model.  Later  we  will  demonstrate  the  technique  using  images  that  follow  a 
negative  binomial  noise  model.  Ultimately,  the  technique  should  work  regardless  of 
the  noise  distribution  assuming  the  correct  iterative  deblurring  algorithm  and  con¬ 
vergence  criteria  are  used. 

This  technique  relies  on  detecting  the  point  where  the  MSE  between  the  collected 
data  and  the  image  estimate  is  lower  than  the  estimated  data  variance.  At  this  point 
we  can  assume  that  any  further  iterations  would  only  serve  to  fit  the  estimates  to  the 
noise.  In  other  words,  for  2-D  images,  iterations  would  cease  when 

^^(d(u,v)  -  1  (u,i>))  <^2^2v(u,v),  (3.7) 

U  V  U  V 

where  V  is  the  actual  image  variance.  Assuming  the  collected  MCFA  follows  Poisson 
statistics,  the  deblurring  algorithm  employed  would  be  the  Richardson-Lucy  algo¬ 
rithm  in  (3.2),  and 

V  (u,  v)  —  d  (u,  v) .  (3.8) 
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The  relationship  in  (3.8)  is  allowed  since  the  assumption  can  be  made  that  the  inten¬ 
sity  captured  by  each  detector  is  independent  of  other  detectors,  and  each  intensity 
can  essentially  be  thought  of  as  an  independent  Poisson  random  variable.  The  sum¬ 
mation  of  these  random  variables  is  therefore  a  good  approximation  for  the  total 
image  variance. 

Due  to  the  photon  counting  nature  of  many  imaging  applications,  the  Poisson  dis¬ 
tribution  is  often  employed  as  a  statistical  distribution  for  the  detected  images.  How¬ 
ever,  due  to  the  highly  coherent  nature  of  laser  light,  images  detected  by  a  LADAR 
sensor  often  follow  more  of  a  negative  binomial  distribution.  Fortunately,  the  robust¬ 
ness  of  the  CoV  technique  allows  for  employment  in  this  scenario  as  well.  MacDonald 
derived  an  iterative  MLE  where  the  noise  is  dominated  by  laser  speckle  as 

U  V  x  7 

U  V  x  7 


(3.9) 


o(x,y)mw  =  o(x,y) 


where  Ai  is  the  coherence  parameter  of  the  light  [44],  Using  the  deblurring  algorithm 
in  (3.9),  we  would  again  iterate  until  the  relationship  in  (3.7)  is  satisfied  where 


V  ( u ,  v)  =  d  ( u ,  v) 


d  (u,  v) 

M 


(3.10) 


The  diagram  in  Figure  3.3  demonstrates  how  this  technique  could  be  employed  in 
an  operational  scenario  where  processing  time  and  computational  power  are  limited. 
In  this  scenario,  images  are  collected  and  fed  into  the  ro  estimation  process.  At  any 
point  in  time,  the  best  possible  estimate  for  r0  can  be  drawn  upon  for  deblurring 
an  image.  However,  in  parallel,  the  r0  estimation  loop  will  continue  to  work  on 
characterizing  the  current  atmospheric  seeing  conditions.  One  of  the  key  advantages 
to  this  algorithm  is  that  it  is  easy  to  parallelize.  Even  with  a  common  home  computer 
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Image  fed  to  Outer  Loop 


Collected  Image 


r0  estimate 
requested  by 
user 


Figure  3.3:  Potential  employment  scenario  for  the  CoV  technique,  where  the  ro  esti¬ 
mation  loop  is  allowed  to  continually  execute  on  recently  collected  images.  The  most 
recent  estimate  for  ro  can  then  be  fed  to  an  iterative  deblurring  algorithm  to  provide 
rapid  results  to  the  user. 


that  has  a  multi-core  processor,  it  is  possible  to  simultaneously  test  multiple  values  of 
ro  for  convergence.  The  r o  to  employ  in  the  deblurring  algorithm  would  be  the  lowest 
value  of  r0  for  which  (3.7)  is  satisfied.  As  MacManus  pointed  out  and  as  indicated 
in  Figure  3.4,  we  can  also  further  enhance  the  process  by  first  conducting  a  coarse 
estimate  of  ro  followed  by  a  more  refined  estimate  [45].  MacDonald  references  the 
employment  of  a  convergence  test  for  ceasing  iterations  in  his  algorithm;  however,  it 
is  apparent  that  he  did  not  utilize  this  test  as  a  sufficient  criteria  for  identifying  the 
correct  value  of  atmospheric  seeing  [41].  As  previously  mentioned,  if  allowed  sufficient 
time,  the  relationship  in  (3.7)  will  be  satisfied  for  the  correct  r0  and  all  values  higher. 
Given  adequate  SNR  in  the  measured  images,  the  criteria  will  never  be  attained  for 
low  estimates  of  r0.  The  following  sections  will  demonstrate  the  employment  of  this 
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Coarse  Search 


lcm  2  cm  3  cm 


4  cm 


— 

4.3  cm 

- / 


Fine  Search 


Figure  3.4:  In  this  demonstration,  the  true  value  for  ro  is  4.3  cm.  However,  the 
algorithm  first  converges  at  5  cm  using  a  1  cm/step  course  search.  It  then  accomplishes 
a  0.1  cm/step  fine  search  to  converge  to  the  true  value  at  4.3  cm. 


technique  on  images  with  Poisson  and  negative  binomial  noise  to  show  that  a  priori 
information  is  not  required  to  achieve  accurate  estimates  for  tq. 


3.2  Simulation 

The  following  results  will  demonstrate  the  utility  of  the  CoV  technique  and  com¬ 
pare  the  results  to  the  MAP  algorithm  developed  by  MacDonald  for  images  with 
Poisson  and  negative  binomial  noise.  The  optical  specifications  listed  in  Table  3.1 
and  used  for  the  simulations  were  not  limited  to  what  could  readily  be  obtained  for 
experimentation.  Rather,  the  specifications  were  chosen  to  mimic  what  could  po¬ 
tentially  be  incorporated  into  a  targeting  pod  design  based  on  size  limitations.  The 
specifications  will  allow  for  properly  sampled  images  according  to  (2.7).  We  will  first 
consider  simulation  results  using  a  fully  illuminated  scene.  We  will  then  simulate 
conditions  for  astrophotography  where  the  scenes  are  only  partially  illuminated. 
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Table  3.1:  Simulated  System  Specifications. 


Parameter  Name 

Defined  Value 

Mean  wavelength  (A) 
Pixel  size 

Sensor  focal  length  (/)) 
Aperture  diameter  ( D ) 
Coherence  Parameter  (iff) 

600  nm 

5  /an  x  5  /an 

3  m 

15  cm 

30 

3.2.1  Fully  Illuminated  Scenes. 

In  this  section  we  will  present  results  from  images  that  are  fully  illuminated.  By 
fully  illuminated,  we  mean  that  the  overwhelming  majority  of  the  scene  is  not  dark, 
and  contains  varying  levels  of  contrast.  Recall  that  MacDonald’s  algorithm  defined  a 
prior  for  r0  of  exponential  form  scaled  by  the  total  number  of  pixels  in  the  detector 
as  shown  in  (3.6). 

In  the  following  examples,  the  cameraman  photo  built  into  MATLAB®  is  blurred 
using  a  total  OTF  that  is  the  product  of  the  diffraction  limited  OTF  and  an  aver¬ 
age  short-exposure  OTF  with  various  levels  of  r'o •  Multiple  trials  will  be  conducted 
with  MCFA  images  composed  of  1,  10,  20,  30,  50  and  100  individual  frames  with 
independent  realizations  of  Poisson  noise  to  demonstrate  the  effects  of  SNR  on  each 
algorithm.  The  original,  blurred  and  an  example  of  a  recovered  image  are  shown  in 
Figure  3.5.  In  order  to  implement  MacDonald’s  algorithm,  we  either  need  an  initial 
estimate  on  the  average  value  for  atmospheric  seeing,  ravg ,  or  we  can  initialize  it  to 
the  aperture  diameter  if  no  estimate  can  be  made.  For  purposes  of  fair  comparison, 
we  will  assume  that  no  prior  estimates  are  known  for  atmospheric  seeing,  and  ravg 
will  be  initialized  to  the  aperture  diameter. 

Table  3.2  shows  that  the  CoV  technique  and  MacDonald’s  algorithm  produce 
nearly  identical  results  with  the  only  exceptions  highlighted  in  bold.  MCFA  images 
consisting  of  more  frames  take  longer  to  converge  due  to  the  higher  intensities  at  each 
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True  Image 


Blurred/Noisy  Image 


Image  Estimate 


Figure  3.5:  In  this  demonstration,  the  true  image  (left)  is  blurred  with  an  average 
short-exposure  OTF  with  an  r0  of  2.6  cm.  The  blurred/noisy  image  (center)  is  the 
summation  of  30  individual  frames  with  independent  realizations  of  Poisson  noise.  The 
image  estimate  (right)  was  obtained  using  the  best  estimate  of  ?’o=2.6  cm  with  the  cap 
on  the  number  of  iterations  set  to  5000  for  the  CoV  technique. 


pixel  associated  with  the  summation  of  individual  frames.  While  the  algorithm  does 
not  always  converge  to  the  true  value  of  r0,  the  estimated  value  was  always  within 
0.5  cm  of  the  true  value.  The  estimates  for  r'o  often  appear  to  be  lower  than  the 
true  value,  and  this  is  likely  due  to  the  algorithm’s  attempt  to  remove  minor  focus 
errors  that  were  not  accounted  for  when  assigning  the  image  as  truth  data.  The 
PSF  arising  from  minor  focus  error  will  blur  the  image  in  a  way  that  may  not  be 
easily  distinguished  from  an  atmospheric  blur  [70].  Therefore  minor  focus  error  could 
contribute  to  the  low  estimates  for  r$.  Additionally,  a  large  portion  of  the  scene 
consists  of  background  objects  that  are  likely  affected  by  a  lower  ro  than  what  is  in 
the  foreground.  Therefore  the  CoV  technique  is  likely  attempting  to  account  for  this. 
This  assessment  is  drawn  from  and  supported  by  the  fact  that,  provided  we  have  an 
adequate  SNR,  the  ideally  simulated  scenes  in  Section  3.2.2  never  converge  below  the 
true  value  regardless  of  the  number  of  iterations  the  algorithm  is  allowed  to  perform. 
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Table  3.2:  Results  for  Fully  Illuminated  Image  Simulation  with  Poisson  Noise  (True 
r o  =  2.6  cm).  Results  in  Bold  Indicate  the  Algorithm  that  Performed  Worse  for  a 
Particular  Scenario 


Max  Iterations  Allowed  -  1000 

Frames 

SNR  (dB) 

MAP  Estimator  (cm) 

CoV  (cm) 

Result 

Error 

Result 

Error 

1 

41.5 

2.4 

-0.2 

2.4 

-0.2 

10 

51.5 

2.7 

0.1 

2.7 

0.1 

20 

54.5 

2.8 

0.2 

2.8 

0.2 

30 

56.2 

2.9 

0.3 

2.9 

0.3 

50 

58.4 

3.0 

0.4 

3.0 

0.4 

100 

61.5 

3.2 

0.6 

3.1 

0.5 

Max  Iterations  Allowed  -  500( 

3 

Frames 

SNR  (dB) 

MAP  Estimator  (cm) 

CoV  (cm) 

Result 

Error 

Result 

Error 

1 

41.5 

2.1 

-0.5 

2.2 

-0.4 

10 

51.5 

2.4 

-0.2 

2.4 

-0.2 

20 

54.5 

2.5 

-0.1 

2.5 

-0.1 

30 

56.2 

2.5 

-0.1 

2.6 

0 

50 

58.4 

2.6 

0 

2.6 

0 

100 

61.5 

2.7 

0.1 

2.7 

0.1 

Max  Iterations  Allowed  -  10000 

Frames 

SNR  (dB) 

MAP  Estimator  (cm) 

CoV  (cm) 

Result 

Error 

Result 

Error 

1 

41.5 

2.1 

-0.5 

2.1 

-0.5 

10 

51.5 

2.3 

-0.3 

2.3 

-0.3 

20 

54.5 

2.4 

-0.2 

2.4 

-0.2 

30 

56.2 

2.4 

-0.2 

2.4 

-0.2 

50 

58.4 

2.5 

-0.1 

2.5 

-0.1 

100 

61.5 

2.6 

0 

2.6 

0 

Max  Iterations  Allowed  -  2000 

i0 

Frames 

SNR  (dB) 

MAP  Estimator  (cm) 

CoV  (cm) 

Result 

Error 

Result 

Error 

1 

41.5 

2.1 

-0.5 

2.1 

-0.5 

10 

51.5 

2.3 

-0.3 

2.3 

-0.3 

20 

54.5 

2.3 

-0.3 

2.3 

-0.3 

30 

56.2 

2.3 

-0.3 

2.4 

-0.2 

50 

58.4 

2.4 

-0.2 

2.4 

-0.2 

100 

61.5 

2.5 

-0.1 

2.5 

-0.1 
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The  performance  of  the  CoV  technique  is  based  on  the  quality  of  the  blurred  image 
and  the  amount  of  time  available  for  processing.  Provided  enough  time  is  allowed, 
(3.7)  will  be  satisfied  for  the  best  estimate  of  tq.  However,  allowing  too  much  time 
does  not  present  a  problem  for  this  algorithm.  By  observing  the  difference  between 
the  left-side  and  right-side  of  (3.7)  at  each  iteration,  we  can  update  an  estimate  for 
the  number  of  remaining  iterations  required  for  convergence,  El,  using 

IV  (n)  —  Y,  (d(u,v)  -  I  (u,v))2 

u,v= 1 

BV=  £  V(u,v)  ,  (3.11) 

u,v=l 

FT—  IV(n—l)—BV 
~  IV(n—l)—IV(n) 

where  n  represents  the  iteration  number,  BV  is  the  variance  of  the  collected  image, 
and  IV  is  the  MSE  between  the  collected  images  and  the  non-noisy  estimate.  Es¬ 
sentially  the  relationships  in  (3.11)  are  used  to  predict  how  long  the  algorithm  will 
have  to  iterate  based  on  the  current  rate  of  convergence  [45].  Figure  3.6  demonstrates 
that  when  the  value  of  ro  is  too  low,  the  estimated  number  of  remaining  iterations 
diverges. 
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Figure  3.6:  Estimated  iterations  remaining  for  an  MCFA  image  composed  of  30  inde¬ 
pendent  frames,  (a)  Coarse  estimation  shows  convergence  for  rg  values  greater  than  3 
cm,  but  divergence  for  values  of  2  cm  or  less  when  the  true  ro—2.6  cm.  (b)  Fine  estima¬ 
tion  with  a  cap  of  5000  iterations  shows  convergence  for  ro  values  of  2.6  cm  or  greater. 
Based  on  experience  with  this  algorithm,  it  is  expected  that  convergence  will  occur  for 
an  ro  value  of  2.4  cm  due  to  the  concave  down  nature  of  the  curve  as  supported  by  the 
results  in  Table  3.2. 
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We  now  repeat  this  experiment  in  the  presence  of  negative  binomial  noise  to  sim¬ 
ulate  the  expected  results  from  laser  illuminated  imagery.  Figure  3.7  again  demon¬ 
strates  the  results  using  an  MCFA  consisting  of  30  frames  and  a  cap  on  the  maximum 
iterations  set  to  5000.  A  close  inspection  of  the  blurred/noisy  image  reveals  a  sig¬ 
nificant  and  visible  increase  in  overall  noise  variance.  As  expected,  this  does  impact 
the  final  results.  Table  3.3  summarizes  the  results  obtained  for  the  CoV  technique, 
as  well  as  the  MAP  estimator  using  the  introduction  of  a  prior  for  the  distribution  of 
ro.  However,  in  the  case  of  negative  binomial  noise,  the  differences  in  the  results  are 
more  significant.  By  introducing  the  prior,  the  tendency  to  underestimate  r0  is  more 
pronounced.  This  trend  of  underestimation  of  tq  was  also  noticed  by  MacDonald  [41] . 
Again,  it  is  expected  that  focus  error  in  the  original  image  is  a  contributing  factor 
in  the  underestimation  of  r0.  This  presents  an  interesting  topic  for  potential  future 
research. 


True  Image 


Blurred/Noisy  Image 


Image  Estimate 


Figure  3.7:  In  this  demonstration,  the  true  image  (left)  is  blurred  with  an  average 
short-exposure  OTF  with  an  ro  of  2.6  cm.  The  blurred/noisy  image  (center)  is  the 
summation  of  30  individual  frames  with  independent  realizations  of  negative  binomial 
noise.  The  image  estimate  (right)  was  obtained  using  the  best  estimate  of  r0=2.6  cm 
with  the  cap  on  the  number  of  iterations  set  to  5000  for  the  CoV  technique. 
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Table  3.3:  Results  for  Fully  Illuminated  Image  Simulation  with  Negative  Binomial 
Noise  (True  rg  =  2.6  cm).  Results  in  Bold  Indicate  the  Algorithm  that  Performed 
Worse  for  a  Particular  Scenario 


Max  Iterations  Allowed  -  1000 

Frames 

SNR  (dB) 

MacDonald’s  Algorithm  (cm) 

CoV  (cm) 

Result 

Error 

Result 

Error 

1 

14.7 

0.9 

-1.7 

1.5 

-1.1 

10 

24.7 

2.1 

-0.5 

2.5 

-0.1 

20 

27.7 

2.4 

-0.2 

2.8 

0.2 

30 

29.4 

2.6 

0 

3.0 

0.4 

50 

31.6 

3.0 

0.4 

3.1 

0.5 

100 

34.7 

3.2 

0.6 

3.2 

0.6 

Max  Iter  at 

;ions  Allowed  -  5000 

Frames 

SNR  (dB) 

MacDonald’s  Algorithm  (cm) 

CoV  (cm) 

Result 

Error 

Result 

Error 

1 

14.7 

0.9 

-1.7 

1.4 

-1.2 

10 

24.7 

1.8 

-0.8 

2.0 

-0.6 

20 

27.7 

2.0 

-0.6 

2.2 

-0.4 

30 

29.4 

2.2 

-0.4 

2.4 

-0.2 

50 

31.6 

2.4 

-0.2 

2.4 

-0.2 

100 

34.7 

2.5 

-0.1 

2.5 

-0.1 

Max  Iterations  Allowed  -  10000 

Frames 

SNR  (dB) 

MacDonald’s  Algorithm  (cm) 

CoV  (cm) 

Result 

Error 

Result 

Error 

1 

14.7 

0.9 

-1.7 

1.3 

-1.3 

10 

24.7 

1.8 

-0.8 

1.9 

-0.7 

20 

27.7 

2.0 

-0.6 

2.1 

-0.5 

30 

29.4 

2.1 

-0.5 

2.2 

-0.4 

50 

31.6 

2.2 

-0.4 

2.3 

-0.3 

100 

34.7 

2.3 

-0.3 

2.4 

-0.2 

Max  Iterations  Allowed  -  20000 

Frames 

SNR  (dB) 

MacDonald’s  Algorithm  (cm) 

CoV  (cm) 

Result 

Error 

Result 

Error 

1 

14.7 

0.8 

-1.8 

1.3 

-1.3 

10 

24.7 

1.8 

-0.8 

1.9 

-0.7 

20 

27.7 

1.9 

-0.7 

2.0 

-0.6 

30 

29.4 

2.0 

-0.6 

2.0 

-0.6 

50 

31.6 

2.1 

-0.5 

2.2 

-0.4 

100 

34.7 

2.2 

-0.4 

2.2 

-0.4 
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At  this  point,  the  functionality  of  the  CoV  technique  has  been  demonstrated 
for  fully  illuminated  scenes.  Further  experimentation  was  conducted  on  simulated 
stellar  targets  to  identify  the  potential  for  measurement  of  atmospheric  seeing  on 
scenes  where  the  majority  of  the  image  consists  only  of  background  illumination  and 
noise.  Since  these  targets  are  fully  simulated,  and  thus  inherently  perfectly  focused, 
underestimation  of  r0  was  not  expected  to  be  a  problem  for  the  CoV  technique. 
Provided  the  algorithm  is  allowed  enough  time  to  iterate,  convergence  to  the  correct 
ro  should  be  achieved.  However,  the  scaling  factor  of  M 2  in  the  prior  (3.6)  was 
expected  to  still  cause  some  bias  in  the  estimates  for  r0  using  MacDonald’s  algorithm. 

3.2.2  Partially  Illuminated  Scenes. 

Space  Situational  Awareness  (SSA)  is  a  key  mission  of  the  United  States  Air  Force 
Space  Command  and  was  a  key  motivator  for  MacManus’  research.  One  aspect  of 
SSA  involves  using  both  telescope  networks  and  radars  to  detect,  identify,  record  and 
track  all  man-made  objects  orbiting  the  Earth.  Knowing  the  exact  locations  of  these 
orbiting  objects  in  space  is  crucial  for  future  space  operation  safety.  Any  debris  that 
results  from  international  space  operations  will  be  an  ongoing  risk  to  US  assets  for 
years  to  come  as  the  orbits  of  the  debris  degrade  toward  Earth.  The  SSA  mission  only 
increases  in  importance  as  additional  high  value  assets  are  placed  in  orbit.  This  is 
merely  a  single  example  and  justification  of  the  importance  for  deblurring  techniques 
applicable  to  astronomical  images  [45]. 

The  following  simulations  will  consider  three  separate  target  configurations.  We 
will  look  at  a  single  point  source  that  could  be  representative  of  a  star,  a  scene  that 
has  multiple  point  sources  arranged  throughout  the  image  and  finally  we  will  look 
at  a  cross  bar  pattern,  as  in  Figure  3.8.  We  will  again  consider  both  Poisson  and 
negative  binomial  assumptions  with  various  SNR  levels. 


Point  Source 


Multiple  Point  Sources 


Cross  Bar  Pattern 


Figure  3.8:  (left)  This  scene  is  representative  of  a  single  star  or  point  source,  (center) 
This  scene  contains  multiple  point  sources  with  varying  intensities  and  spacings.  The 
spacing  between  the  two  point  sources  in  the  center  of  the  image  is  a  single  pixel, 
(right)  This  scene  contains  a  cross  bar  pattern. 


The  testing  in  Section  3.2.1  revealed  that  both  algorithms  are  impacted  by  SNR 
and  focus  error.  However,  the  simulations  also  revealed  that  no  gain  in  performance 
was  realized  through  the  introduction  of  the  prior  for  to,  and  that  the  CoV  technique 
exhibited  promise  for  estimation  of  tq  as  well  as  a  deblurred  scene.  The  results  that 
follow  add  support  for  the  hypothesis  that  underestimation  of  r0  was  a  function  of 
both  SNR  and  focus  error.  The  SNR  for  the  various  MCFAs  used  in  the  following 
simulations  is  identified  in  Table  3.4.  As  expected,  the  SNR  for  MCFAs  with  Poisson 
noise  is  higher  than  the  SNR  for  MCFAs  with  negative  binomial  noise. 

For  demonstration  purposes,  the  algorithm  was  allowed  a  cap  of  1,000,000  iter¬ 
ations.  Even  under  these  conditions,  underestimation  was  never  a  factor  for  images 
with  adequate  SNR  using  the  CoV  technique.  At  this  point,  it  is  unknown  if  it  would 
be  possible  to  predict  the  precise  level  of  SNR  at  which  the  algorithm  will  converge 
to  the  correct  value  of  r0  since  the  relationship  appears  to  be  scene  or  contrast  depen¬ 
dent.  However,  we  can  conclude  that  higher  levels  for  SNR  will  yield  better  results. 
Additionally,  even  at  low  values  of  SNR  we  achieve  reasonable  estimations  for  the 
deblurred  scene  and  r0  using  the  CoV  technique.  On  the  other  hand,  by  introducing 
a  prior  for  the  distribution  of  ro,  underestimation  is  more  prevalent. 
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Table  3.4:  Signal  to  Noise  Ratio  for  Simulation  Data 


Signal  to  Noise  Ratio  (dB)  for  Poisson  MCFAs 

Frames  in  MCFA 

Point  Source 

Multiple  Point  Sources 

Cross  Bar 

1 

15.0 

17.2 

30.0 

10 

24.8 

26.2 

39.3 

20 

28.0 

29.2 

42.9 

30 

29.6 

30.7 

44.7 

50 

31.9 

33.3 

46.9 

Signal  to  Noise  Ratio  (dB)  for  Negative  Binomial  MCFAs 

Frames  in  MCFA 

Point  Source 

Multiple  Point  Sources 

Cross  Bar 

1 

11.8 

13.1 

14.4 

10 

22.6 

23.7 

25.7 

20 

26.3 

27.0 

29.3 

30 

27.5 

28.1 

30.5 

50 

29.4 

30.0 

32.7 

Tables  3.5  and  3.7  demonstrate  the  utility  of  the  CoV  algorithm  on  partially 
illuminated  scenes  with  Poisson  and  negative  binomial  noise  respectively.  Trials  where 
the  estimate  for  r0  matched  the  true  value  are  shown  in  bold.  From  these  results,  we 
conclude  that  provided  enough  frames  are  properly  registered  and  averaged  to  provide 
adequate  SNR,  and  enough  time  is  allowed  for  convergence,  the  value  of  r0  can  be 
estimated  to  within  1  mm  for  the  simulated  optical  configuration.  In  cases  where  we 
have  low  SNR  the  algorithm  will  tend  to  underestimate,  but  this  is  expected  since  the 
noise  power  in  the  images  masks  some  of  the  high  frequency  content.  If  insufficient 
time  is  allowed  for  convergence,  the  algorithm  will  produce  a  high  estimate  for  ro,  as 
observed  when  the  cap  for  iterations  was  limited  to  1,000.  Additionally,  we  notice 
that  in  cases  of  adequate  SNR,  we  do  not  have  a  problem  of  underestimation  of  r0 
since  focus  error  was  not  present  in  these  images. 
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Table  3.5:  Results  for  Partially  Illuminated  Image  Simulation  with  Poisson  Noise  (True 
r o  =  3.3  cm).  Results  in  Bold  Indicate  the  Estimate  Matched  the  True  Value  for  rp. 


Max  Iterations  Allowed  -  1000 

Frames 

MAP  Estimator  (cm) 

CoV  (cm) 

Point 

Multi-Point 

Cross  Bar 

Point 

Multi-Point 

Cross  Bar 

1 

0.1 

1.5 

2.6 

3.2 

3.2 

3.4 

10 

2.9 

3.1 

3.3 

3.3 

3.5 

4.2 

20 

3.2 

3.2 

3.4 

3.3 

3.6 

4.6 

30 

3.3 

3.3 

3.5 

3.3 

3.8 

5.1 

50 

3.3 

3.3 

3.6 

3.4 

3.9 

5.9 

Max  Iterations  Allowed  -  1C 

o 

o 

o 

Frames 

MAP  Estimator  (cm) 

CoV  (cm) 

Point 

Multi-Point 

Cross  Bar 

Point 

Multi-Point 

Cross  Bar 

1 

0.1 

1.3 

2.6 

3.2 

3.2 

3.2 

10 

2.9 

3.0 

3.2 

3.3 

3.3 

3.3 

20 

3.1 

3.2 

3.3 

3.3 

3.3 

3.3 

30 

3.2 

3.2 

3.3 

3.3 

3.3 

3.3 

50 

3.2 

3.2 

3.3 

3.3 

3.3 

3.4 

Max  Iterations  Allowed  -  20000 

Frames 

MAP  Estimator  (cm) 

CoV  (cm) 

Point 

Multi-Point 

Cross  Bar 

Point 

Multi-Point 

Cross  Bar 

1 

0.1 

1.3 

2.6 

3.2 

3.2 

3.2 

10 

2.9 

3.0 

3.2 

3.3 

3.3 

3.3 

20 

3.1 

3.1 

3.3 

3.3 

3.3 

3.3 

30 

3.1 

3.2 

3.3 

3.3 

3.3 

3.3 

50 

3.2 

3.2 

3.3 

3.3 

3.3 

3.3 

Max  Iterations  Allowed  -  1000000 

Frames 

MAP  Estimator  (cm) 

CoV  (cm) 

Point 

Multi-Point 

Cross  Bar 

Point 

Multi-Point 

Cross  Bar 

1 

0.1 

1.3 

2.5 

3.2 

3.2 

3.2 

10 

2.8 

3.0 

3.2 

3.3 

3.3 

3.3 

20 

3.1 

3.1 

3.3 

3.3 

3.3 

3.3 

30 

3.1 

3.2 

3.3 

3.3 

3.3 

3.3 

50 

3.2 

3.2 

3.3 

3.3 

3.3 

3.3 
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In  Figure  3.9  we  demonstrate  the  performance  of  the  CoV  technique  for  the  three 
target  types  in  the  presence  of  Poisson  noise.  Using  the  RL  deconvolution  algorithm 
with  a  cap  on  the  number  of  iterations  set  to  10,000,  and  an  MCFA  consisting  of 
30  frames,  we  were  able  to  recover  the  correct  r0  in  all  cases.  For  the  point  source 
targets,  the  blurring  effects  of  the  simulated  atmosphere  reduce  the  intensity  to  the 
point  that  it  is  difficult  to  visually  identify  the  various  point  sources.  However,  when 
deconvolution  is  completed,  each  of  the  sources  is  easily  identified.  For  the  cross  bar 
pattern,  the  process  of  deconvolution  makes  it  much  easier  to  identify  the  structure 
of  the  target  pattern.  While  it  may  seem  that  a  cap  of  10,000  iterations  is  unrea¬ 
sonable,  the  algorithm  will  only  take  as  much  time  as  needed  to  converge  within  this 
upper  bound.  For  instance,  if  the  termination  criteria  is  met  before  the  upper  bound 
on  iterations  is  achieved,  the  algorithm  will  terminate.  Even  with  a  standard  home 
computer  with  a  2.7GHz  Intel®  Core™  i5  processor  and  16GB  of  memory,  conver¬ 
gence  was  achieved  in  a  reasonably  short  period  of  time  for  all  target  types  as  shown 
in  Table  3.6.  In  this  example,  coarse  convergence  is  to  the  nearest  centimeter  and 
Table  3.6:  Convergence  Times  for  Simulations  Shown  in  Figure  3.9. 


Target  Type 

Coarse  Convergence 

Fine  Convergence 

Total 

Iterations 

Time  (sec) 

Iterations 

Time  (sec) 

Time  (sec) 

Single  Point 

627 

1.66 

968 

2.49 

4.15 

Multi-Point 

800 

1.90 

1,709 

4.45 

6.35 

Cross  Bar 

2,123 

4.26 

8,275 

17.51 

21.77 

fine  convergence  goes  to  the  nearest  millimeter.  The  algorithm  could  be  manually 
interrupted  sooner  if  needed.  At  which  point,  the  lowest  value  of  r0  that  has  allowed 
convergence  would  be  returned  as  the  answer.  For  instance,  in  the  example  shown  in 
Figure  3.9,  if  we  interrupted  the  routine  after  4.26  seconds  for  the  cross  bar  pattern, 
we  would  get  an  estimated  r0  of  4  cm.  The  estimate  continues  to  be  refined  until  the 
best  estimate  is  achieved  after  21.77  seconds. 
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True  Image 


Blurred/Noisy  Image 


Image  Estimate 


(a)  Single  Point  Source 


(b)  Multiple  Point  Sources 


Image  Estimate 


True  Image 


Blurred/Noisy  Image 


Image  Estimate 


(c)  Cross  Bar  Pattern 


Figure  3.9:  Demonstration  of  the  CoV  technique  using  the  RL  deconvolution  algorithm 
with  an  MCFA  consisting  of  30  frames  for  the  single  point  source  (a),  multiple  point 
sources  (b),  and  the  cross  bar  pattern  (c). 
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Table  3.7:  Results  for  Partially  Illuminated  Image  Simulation  with  Negative  Binomial 
Noise  (True  ro  =  3.3  cm).  Results  in  Bold  Indicate  the  Estimate  Matched  the  True 
Value  for  r(J . 


Max  Iterations  Allowed  -  1000 

Frames 

MAP  Estimator  (cm) 

CoV  (cm) 

Point 

Multi-Point 

Cross  Bar 

Point 

Multi-Point 

Cross  Bar 

1 

0.1 

1.0 

1.8 

3.1 

3.0 

5.0 

10 

2.4 

2.7 

2.8 

3.4 

3.7 

5.1 

20 

2.9 

3.0 

3.0 

3.5 

4.0 

5.8 

30 

3.0 

3.1 

3.1 

3.7 

4.8 

6.5 

50 

3.2 

3.2 

3.2 

3.9 

5.0 

8.5 

Max  Iterations  Allowed  -  10000 

Frames 

MAP  Estimator  (cm) 

CoV  (cm) 

Point 

Multi-Point 

Cross  Bar 

Point 

Multi-Point 

Cross  Bar 

1 

0.1 

0.5 

1.7 

3.1 

2.9 

3.4 

10 

2.3 

2.6 

2.7 

3.3 

3.2 

3.4 

20 

2.7 

2.9 

3.0 

3.3 

3.3 

3.5 

30 

2.9 

3.0 

3.1 

3.3 

3.3 

3.7 

50 

3.0 

3.1 

3.3 

3.3 

3.3 

4.0 

Max  Iterations  Allowed  -  2C 

o 

o 

o 

Frames 

MAP  Estimator  (cm) 

CoV  (cm) 

Point 

Multi-Point 

Cross  Bar 

Point 

Multi-Point 

Cross  Bar 

1 

0.1 

0.4 

1.7 

3.1 

2.9 

3.2 

10 

2.3 

2.6 

2.7 

3.2 

3.2 

3.2 

20 

2.7 

2.9 

3.1 

3.2 

3.3 

3.3 

30 

2.9 

3.0 

3.1 

3.3 

3.3 

3.3 

50 

3.0 

3.1 

3.2 

3.3 

3.3 

3.3 

Max  Iterations  Allowed  -  IOC 

o 

o 

o 

o 

Frames 

MAP  Estimator  (cm) 

CoV  (cm) 

Point 

Multi-Point 

Cross  Bar 

Point 

Multi-Point 

Cross  Bar 

1 

0.1 

0.4 

1.7 

3.1 

2.9 

3.2 

10 

2.3 

2.5 

2.7 

3.2 

3.2 

3.2 

20 

2.7 

2.8 

3.0 

3.2 

3.3 

3.3 

30 

2.8 

3.0 

3.1 

3.3 

3.3 

3.3 

50 

3.0 

3.1 

3.1 

3.3 

3.3 

3.3 
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3.3  Experimental  Results 


The  optical  configuration  shown  in  Figure  2.3  was  used  to  obtain  properly  sampled 
images.  The  specifications  for  this  setup  are  listed  in  Table  2.1.  The  experiments  for 
fully  illuminated  scenes  will  use  the  image  in  Figure  3.10  as  a  target.  This  target  will 
be  placed  indoors,  10  meters  from  the  sensor  where  turbulence  and  illumination  can  be 
controlled.  The  incorporation  of  a  step  in  the  bottom  portion  of  the  scene  will  permit 
the  measurement  of  the  true  r0  for  comparison  with  the  estimated  values.  In  order  to 
create  a  turbulent  atmosphere  to  image  through,  a  heat  source  was  directed  in  front 
of  the  telescope  aperture.  This  technique  allowed  for  the  generation  of  repeatable 
turbulent  atmospheres  with  r0  values  in  the  sub-centimeter  range.  Without  the  use 
of  this  heat  source,  all  of  the  images  would  likely  have  been  at  or  near  the  diffraction 
limit  making  validation  of  the  CoV  technique  difficult. 


E  E  E  E 
Q  Q  qq 
DDdd 
P  P  p  p 
C  CCc 


Figure  3.10:  Scene  used  for  each  of  the  fully  illuminated  experiments.  The  top  portion 
of  the  scene  includes  various  characters  of  decreasing  size.  The  bottom  portion  of  the 
scene  has  a  step  target  to  allow  for  measurement  of  the  true  tq. 


3.3.1  Fully  Illuminated  Scenes  -  Natural  Light. 

In  the  following  two  experiments,  the  remote  scene  was  fully  illuminated  by  nat¬ 
ural  lighting.  Given  that  the  light  source  was  incoherent,  the  Poisson  statistical 
distribution  for  the  photon  arrival  rate  applies  [21],  We  will  use  the  Richardson- 
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Lucy  deconvolution  algorithm  (3.2)  with  PSFs  parameterized  by  a  range  of  r0  values 
from  0.1  cm  up  to  the  aperture  diameter  of  5  cm.  The  images  were  blurred  using  a 
heat  source  to  create  various  levels  of  atmospheric  turbulence  and  subsequent  image 
blurring. 

In  Figure  3.11  we  demonstrate  the  ability  to  measure  ro  by  measuring  the  step 
response  from  the  collected  MCFA  as  described  in  Section  2.5.  The  impulse  response 
is  then  found  by  taking  the  derivative  of  this  measured  step  response.  Once  we 
have  the  impulse  response,  we  can  vary  r o  per  the  relationship  in  (2.9)  to  find  the 
theoretical  total  PSF,  htot,  that  minimizes  the  error  between  the  measured  impulse 
response  and  the  theoretical  impulse  response. 


76 


(a)  Collected  MCFA 


(b)  Mean  Step  Response 


(c)  Theoretical/Experimental  OTF  Comparison 

Figure  3.11:  Using  the  step  in  the  bottom  portion  of  the  colected  remote  scene  (a) 
we  can  compute  the  mean  step  response  for  the  collected  MCFA  (b).  From  this  step 
response  we  compute  the  experimental  OTF  and  find  the  theoretical  short  exposure 
OTF  that  minimizes  MSE  between  the  two  (c). 
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In  the  first  experiment  shown  in  Figure  3.12,  r0  was  measured  at  0.4  cm.  Using 
the  CoV  technique,  we  obtain  an  estimate  of  r'o  =  0.5  cm  for  an  error  in  estimation 
of  only  1  mm.  At  this  point  it  is  important  to  note  that  the  edges  of  the  deblurred 
images  are  distorted  due  to  the  implementation  of  the  RL  algorithm  using  the  2-D 
FFT  in  MATLAB.  This  implementation  was  chosen  in  order  to  decrease  the  time 
required  for  execution  and  results  in  a  fair  substitution  as  long  as  the  image  edges 
are  not  of  importance.  Therefore,  when  computing  the  variance  per  the  relationship 
in  (3.7),  the  variance  for  the  image  edges  were  ignored.  In  the  collected  MCFA,  it  is 
difficult  to  discern  the  smaller  font  sizes.  However,  when  deconvolution  is  conducted 
with  an  OTF  parameterized  by  the  estimate  for  ro,  it  is  possible  to  identify  each  of 
the  characters.  In  the  second  experiment  shown  in  Figure  3.13,  r0  was  measured  at 
1.1  cm.  Using  the  CoV  technique,  we  obtain  an  estimate  of  r o  =  1.1  cm.  While 
the  blurring  due  to  turbulence  was  less  severe  in  this  experiment,  improvement  in 
sharpness  of  the  characters  is  again  noted  when  deconvolution  was  conducted  using 
the  estimate  for  ro- 


(a)  Collected  MCFA  (b)  Deblurred  Scene 

Figure  3.12:  (a)  Collected  MCFA  consisting  of  100  registered  frames,  each  with  an 

exposure  time  of  0.001s.  (b)  Deblurred  image  using  the  lowest  ?’o  for  which  convergence 
was  achieved  (ro  =  0.5  cm). 
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(a)  Collected  MCFA  (b)  Deblurred  Scene 

Figure  3.13:  (a)  Collected  MCFA  consisting  of  100  registered  frames,  each  with  an 

exposure  time  of  0.001s.  (b)  Deblurred  image  using  the  lowest  ?’o  for  which  convergence 
was  achieved  (ro  =  1.1  cm). 

3.3.2  Fully  Illuminated  Scenes  -  Laser  Illumination. 

In  the  following  two  experiments,  the  remote  scene  was  illuminated  using  a  laser 
with  a  wavelength  of  630  nm.  Given  that  the  light  source  was  partially  coherent,  with 
a  measured  coherence  parameter,  A4  =  10,  the  negative  binomial  model  for  speckle 
noise  applies  [21].  We  will  use  the  ML  estimator  in  (3.9)  with  PSFs  parameterized 
by  a  range  of  ro  values  from  0.1  cm  up  to  the  aperture  diameter  of  5  cm. 

In  the  first  experiment  shown  in  Figure  3.14,  ro  was  measured  at  0.5  cm.  Using 
the  CoV  technique,  we  obtain  an  estimate  of  r o  =  0.5  cm.  In  the  collected  MCFA, 
it  is  difficult  to  discern  the  smaller  font  sizes,  and  based  on  where  the  illumination 
spot  was  centered,  the  top  two  rows  of  text  are  nearly  illegible.  However,  when 
deconvolution  is  conducted  with  an  OTF  parameterized  by  the  estimate  for  ro  it  is 
possible  to  identify  most  of  the  characters.  It  is  much  easier  to  identify  the  row  of  Qs, 
and  the  top  row  of  Es  is  faintly  visible.  In  the  second  experiment  shown  in  Figure 
3.15,  ro  was  measured  at  1.1  cm.  Using  the  CoV  technique,  we  obtain  an  estimate  of 
rQ  =  1.1  cm.  While  the  blurring  due  to  turbulence  was  less  severe  in  this  experiment, 
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significant  improvement  is  again  noted  when  deconvolution  was  conducted  using  the 
estimate  for  r q. 


Table  3.8:  Summary  of  Results  for  Fully  Illuminated  Image  Experiments. 


Trial 

Estimated  ro  (cm) 

Measured  ro  (cm) 

Natural  Light  -  Low  r0 

0.5 

0.4 

Natural  Light  -  High  r'o 

1.1 

1.1 

Laser  Illumination  -  Low  ro 

0.5 

0.5 

Laser  Illumination  -  High  r0 

1.1 

1.1 

80 


Figure  3.14:  (a)  Collected  MCFA  consisting  of  100  registered  frames,  each  with  an 

exposure  time  of  0.005s.  (b)  Deblurred  image  using  the  lowest  r0  for  which  convergence 
was  achieved  (ro  =  0.5  cm). 


Figure  3.15:  (a)  Collected  MCFA  consisting  of  100  registered  frames,  each  with  an 

exposure  time  of  0.005s.  (b)  Deblurred  image  using  the  lowest  ro  for  which  convergence 
was  achieved  (ro  =  1.1  cm). 
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3.3.3  Correlation  Technique  for  Measurement  of  Atmospheric  Seeing. 


With  the  previous  experiments,  the  true  value  of  r0  could  be  measured  by  imaging 
a  step  target  and  taking  the  derivative  to  find  the  impulse  response  as  shown  in  Fig¬ 
ure  3.11.  The  step  target  could  be  placed  in  line  with  the  desired  remote  scene  such 
that  the  measurements  were  made  through  nearly  identical  columns  of  turbulent  air. 
However,  when  trying  to  measure  the  true  r'o  of  stellar  images  for  comparison  with 
the  experimental  results,  this  was  not  a  viable  solution.  We  could  average  many  short 
exposure  images  of  a  single  star  in  order  to  get  the  average  short  exposure  OTF  pa¬ 
rameterized  by  r0,  however  this  requires  precise  tilt  removal  and  any  shift  estimation 
errors  will  appear  as  attenuation  of  the  short  exposure  OTF  and  underestimation  of 
ro.  Trying  to  average  enough  frames  to  achieve  the  long  exposure  OTF  in  order  to 
find  r0  would  likely  take  thousands  of  images  to  converge  upon  the  optimal  solution. 
Based  on  the  frame  rate  of  the  experimental  collection  system,  it  is  possible  that  the 
value  for  ro  would  change  in  the  time  required  to  gather  this  amount  of  data.  There¬ 
fore,  we  considered  an  alternative  technique  for  measuring  the  value  of  atmospheric 
seeing  that  considers  the  cross  correlation  of  the  collected  images. 

By  considering  all  possible  combinations  of  the  cross  correlations  between  a  series 
of  individual  short  exposure  images,  {Sk  :  k  =  1, ...,  K},  taken  of  a  star,  we  can  find 
the  Cross  Power  Spectral  Density  (CPSD),  Ps(z/X,  z/y),  where 

Ps  (ux,  uy)  =  E  (x,y)}F  {Sk/  (x,y)}]Vk  ^  k' .  (3.12) 

In  other  words,  the  CPSD  is  the  correlation  between  the  normalized  Fourier  trans¬ 
forms  of  all  possible  image  combinations  [30].  For  a  sequence  of  K  images,  there  are 
a  total  of  K (K  —  l)/2  non  redundant  cross  correlations.  Additionally,  the  CPSD  of 
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the  blurred  point  source  can  be  shown  to  have  the  following  relationship 

ps{vx,Vy)  =  E{\H(vx,vy)^}^  (3.13) 

where  E  {\H  (ux,  vy)\2}  is  the  speckle  transfer  function  [59].  Fortunately,  the  speckle 
transfer  function  can  also  be  parameterized  by  r0.  Therefore,  in  order  to  obtain  the 
true  value  for  atmospheric  seeing,  we  can  find  the  value  of  ro  that  minimizes  the  error 
in  (3.13).  At  this  point  we  must  keep  in  mind  that  this  technique  for  measuring  r0  is 
limited  to  cases  where  we  are  imaging  a  point  source.  We  will  use  this  technique  to 
demonstrate  that  the  CoV  technique  will  in  fact  identify  the  correct  ro- 

3.3.4  Partially  Illuminated  Scenes. 

For  the  experiments  involving  partially  illuminated  scenes,  we  chose  to  image  a 
star.  Stars  can  essentially  be  considered  point  sources  of  light,  allowing  us  to  use 
the  technique  presented  in  Section  3.3.3  to  obtain  truth  data  for  comparison  with 
the  CoV  results.  While  the  resultant  deblurred  image  is  intuitive,  the  estimated 
values  for  r0  prove  that  the  technique  can  be  successfully  used  to  measure  r0  for 
partially  illuminated  scenes.  For  the  experiments,  we  chose  a  relatively  bright  star 
near  Polaris  to  image.  This  minimized  the  relative  motion  between  the  FOV  and  the 
imaged  portion  of  the  sky.  All  individual  images  were  taken  using  an  exposure  time  of 
0.001s  to  ensure  that  the  short-exposure  model  was  applicable  [21].  Additionally,  each 
of  the  MCFAs  is  a  compilation  of  20  individual  frames.  Some  experimentation  with 
MCFAs  consisting  of  more  than  20  frames  was  accomplished.  However,  no  increase 
in  performance  was  observed. 

In  Table  3.9,  the  estimated  value  for  ro  was  always  within  0.2  cm  of  the  measured 
value.  Initially  the  cap  on  the  number  of  iterations  was  set  to  1,000.  In  all  cases,  if 
convergence  was  achieved  for  a  particular  level  of  r0,  it  occurred  within  the  first  100 
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Table  3.9:  Results  for  Partially  Illuminated  Image  Experiments. 


Trial 

Estimated  ro  (cm) 

Measured  ro  (cm) 

1 

2.3 

2.1 

2 

2.4 

2.2 

3 

2.3 

2.2 

4 

2.4 

2.1 

5 

2.4 

2.1 

iterations.  We  then  increased  the  number  of  iterations  to  10,000  and  subsequently 
20,000  to  see  if  convergence  could  be  achieved  for  a  lower  level  of  r0,  but  the  results 
remained  the  same  with  the  increased  cap  on  iterations. 


Collected  MCFA  Image  Estimate 


Figure  3.16:  The  collected  MCFA  consisting  of  20  registered  short  exposure  images 
(left)  and  the  deblurred  image  estimate  when  the  average  short  exposure  OTF  with  an 
ro  of  2.4  cm  is  used  for  deconvolution  (right). 

For  the  image  shown  in  Figure  3.16  we  first  conducted  a  course  search  with  a  step 
size  of  1  cm  revealing  that  the  minimum  value  of  r0  that  will  allow  convergence  is  3 
cm.  We  then  step  through  with  a  smaller  step  size  and  show  that  convergence  can 
be  achieved  for  values  of  r0  as  low  as  2.4  cm.  For  this  image,  convergence  occurs 
after  just  68  iterations  for  an  r0  of  2.4  cm.  However,  it  is  divergent  for  anything  less 
than  this  value.  In  Figure  3.16  we  show  the  resultant  deblurred  image  when  using 
this  best  estimate  for  r0  in  conjunction  with  the  RL  deconvolution  algorithm  and  the 
average  short  exposure  OTF.  As  expected  the  image  estimate  approaches  more  of  a 
point  source  than  the  original  image. 


84 


3.4  Extension  of  CoV  Technique  for  3-D  FLASH  LADAR  Images 

Through  the  course  of  this  research,  it  was  found  that  the  accuracy  of  the  CoV 
technique  was  primarily  dependent  on  the  variance  determined  threshold  for  stopping 
the  algorithm.  This  threshold  is  chosen  based  on  the  assumed  noise  model.  For 
FLASH  LADAR  imagery,  this  does  present  a  challenge  since  there  are  again  various 
forms  of  noise  present  in  the  image.  In  order  to  effectively  use  the  CoV  technique 
we  must  be  able  to  either  predict  the  variance  attributed  to  each  source  of  noise,  or 
measure  the  total  variance  through  a  series  of  independent  images. 

Three  dimensional  LADAR  technology  is  receiving  an  increased  interest  as  the 
technology  improves.  Currently,  the  commercially  available  sensors  are  severely 
under-sampled,  and  do  not  experience  the  effects  of  diffraction  from  atmospheric 
turbulence.  However,  as  the  technology  continues  to  progress,  it  is  expected  that 
minimizing  the  effects  of  atmospheric  turbulence  will  be  important.  Conversely,  cer¬ 
tain  applications  such  as  the  imaging  and  tracking  of  space  debris  may  require  an 
optics  configuration  where  the  current  sensors  would  receive  properly  sampled  im¬ 
ages.  In  those  cases,  the  CoV  technique  could  potentially  be  applied  to  identify  the 
PSF  parameterized  by  r0  that  will  deblur  the  scene. 

A  typical  full-waveform  3-D  LADAR  image  is  comprised  of  multiple  2-D  images 
or  frames  separated  by  a  small  time  delta.  Therefore,  the  3-D  image  can  be  flattened 
into  a  2-D  image  by  simply  removing  the  range  information  and  accumulating  the  in¬ 
tensity  information  for  each  pixel  for  the  series  of  individual  frames.  Flattening  a  3-D 
image  into  a  2-D  image  presents  several  challenges  to  the  CoV  technique.  The  pri¬ 
mary  issue  being  that  the  summation  of  multiple  negative  binomial  random  variables 
with  different  means  is  not  another  negative  binomial  random  variable.  Therefore 
it  would  be  difficult  to  determine  the  best  deconvolution  algorithm  to  use  based  on 
the  distribution  of  the  data.  By  applying  the  Central  Limit  Theorem,  it  may  be 
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possible  to  make  a  case  for  employing  a  deconvolution  algorithm  based  on  the  Gaus¬ 
sian  distribution.  Likely  less  of  a  limiting  factor  would  be  the  case  where  the  image 
has  uniform  reflectance  across  a  diversity  of  ranges.  In  this  scenario,  any  contrast 
in  intensity  through  multiple  regions  of  the  image  may  be  lost,  making  it  difficult  to 
determine  the  best  tq.  Fortunately,  the  atmosphere  can  be  considered  static  for  the 
laser  illuminator  pulse  duration  and  subsequent  detector  integration  times  that  are 
common  to  3-D  LADAR  sensors  [21],  [58].  Given  these  challenges,  likely  the  best 
option  for  employing  the  CoV  technique  on  3-D  FLASH  LADAR  images  would  be 
to  treat  each  of  the  2-D  frames  separately.  Based  on  this  premise,  this  technique 
was  originally  explored  as  a  means  of  identifying  the  best  PSF  parameterized  by  ro 
to  be  employed  in  algorithms  such  as  the  multiple  surface  FLASH  LADAR  ranging 
algorithm  that  will  be  discussed  in  Chapter  IV  and  originally  presented  in  [53]. 

Given  the  3-D  FLASH  LADAR  image  model  provided  in  Chapter  II,  we  expect 
the  total  variance  to  consist  of  multiple  components.  Likely  the  most  dominant 
form  of  image  variance  in  illuminated  portions  of  the  scene  will  be  speckle  noise. 
However,  the  additional  detector  noise  or  background  lighting  noise  may  significantly 
impact  the  overall  stopping  threshold,  especially  in  low  illumination  images.  Due  to 
the  independence  of  the  different  noise  sources,  the  total  image  variance  for  the  kth 
frame,  a^k,  will  equal 

°Tk  =  ( aSk  +  aBk)i  (3-14) 

where  a2Sk  is  the  noise  due  to  speckle  illumination  and  a2Bk  is  the  noise  due  to  detector 
bias  for  each  of  the  K  2-D  frames  in  the  3-D  image.  Since  the  3-D  image  is  made 
up  of  numerous  2-D  images  corresponding  with  unique  ranges,  we  need  to  choose  the 
image  that  allows  the  CoV  technique  the  greatest  chance  for  success.  Given  that  the 
target  profile  in  Figure  3.17  (a)  is  centered  within  the  range  gate  of  the  3-D  FLASH 
LADAR  sensor,  we  would  expect  the  individual  2-D  frames  to  have  various  levels  of 


contrast.  Additionally,  it  is  possible  that  the  2-D  frames  corresponding  to  specific 
ranges  may  not  have  a  visible  return  present  as  shown  in  Figure  3.18.  As  shown 
in  Figure  3.17(b),  the  estimates  of  ro  are  generally  within  3  mm  of  the  true  value 
when  the  frames  of  higher  contrast  (frames  6-10  and  13-15)  are  used.  However,  when 
frames  with  low  contrast  are  used,  the  accuracy  of  the  estimates  falls  off.  Additional 
demonstrations  with  the  CoV  technique  on  actual  3-D  FLASH  LADAR  images  will 
be  provided  in  Chapter  V. 
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Figure  3.17:  (a)  Simulated  target  with  a  single  raised  bar  placed  at  a  range  of  103 

meters  and  the  background  is  located  at  105  meters.  The  reflectivity  of  the  single  bar 
is  0.7  and  the  reflectivity  of  the  background  is  1.0.  (b)  Estimated  values  of  ro  for  each 
of  the  2-D  frames  in  the  simulated  3-D  FLASH  LADAR  image  where  the  true  value 
of  r0  =  1.5  cm.  Error  bars  indicate  the  standard  deviation  of  estimates  for  30  separate 
trials. 


Figure  3.18:  Individual  2-D  frames  based  on  the  single  bar  target  in  Figure  3.17  being 
centered  in  the  range  gate  for  a  LADAR  system  with  parameters  listed  in  Table  4.1. 
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3.5  Chapter  Summary 


The  original  focus  of  this  work  was  to  develop  a  blind  deconvolution  technique 
that  could  be  employed  in  a  tactical  military  environment  where  both  time  and  com¬ 
putational  power  are  limited.  The  intent  behind  its  expansion  and  inclusion  in  this 
dissertation  was  to  provide  a  means  of  comparison  with  the  techniques  that  will  be 
provided  in  Chapter  V.  The  CoV  technique  detailed  above  allows  for  rapid  and  ac¬ 
curate  estimations  of  the  atmospheric  OTF  parameterized  by  the  seeing  parameter, 
ro.  As  shown  in  Figure  3.3,  the  technique  can  be  interrupted  after  any  amount  of 
time,  at  which  point  the  best  available  results  would  be  provided.  If  more  time  is  pro¬ 
vided,  the  results  are  generally  enhanced.  Additionally,  the  CoV  technique  reduces 
the  possibility  of  noise  amplification  common  with  iterative  deconvolution  algorithms 
by  ceasing  iteration  once  the  statistically  predicted  variance  is  achieved. 

An  interesting  discovery  was  also  made  through  the  course  of  this  effort  and  is 
highlighted  in  Section  3.2.1.  This  technique  may  be  useful  in  recovering  from  minor 
focus  error  in  the  collected  images.  There  are  similarities  between  the  atmospheric 
OTF  and  the  OTF  that  arises  from  a  minor  focus  error.  As  a  result,  the  estimates 
for  r0  may  be  lower  than  the  true  value.  Therefore,  if  this  algorithm  is  to  be  used  for 
atmospheric  seeing  measurement,  the  images  must  be  in  focus.  Another  source  of  error 
in  estimation  may  arise  if  various  portions  of  the  scene  have  differing  levels  of  r$.  In  the 
cameraman  photo  used  for  the  simulations,  the  objects  in  the  background  were  likely 
impacted  by  a  different  level  of  atmospheric  turbulence  than  those  in  the  foreground. 
A  potential  topic  for  future  research  would  be  the  relationship  between  the  seeing 
parameter  and  focus  interaction  for  varying  levels  of  focus  error  and  atmospheric 
turbulence. 
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IV.  Multiple  Surface  Estimation 


This  chapter  focuses  on  the  development  of  a  multiple  surface  ranging  algorithm 
for  full  waveform  3-D  FLASH  LADAR  images.  In  addition  to  providing  pulse  infor¬ 
mation  for  multiple  surfaces  per  image  pixel,  the  presented  algorithm  also  reduces 
the  effects  of  diffraction  [52],  Simulation  results  will  be  presented  to  demonstrate  the 
utility  of  this  algorithm  in  cases  where  the  blurring  function  is  known.  Experimental 
results  will  be  presented  in  Chapter  V  in  conjunction  with  results  demonstrating  the 
ability  to  solve  for  a  parameterized  blurring  function. 

Currently  there  are  numerous  efficient  techniques  such  as  peak  estimation,  cross¬ 
correlation  and  leading  edge  detection  for  identifying  a  single  surface  per  pixel  [58]. 
Of  these  techniques,  the  cross-correlation  technique  can  be  shown  to  be  extremely  ac¬ 
curate  in  identifying  the  range  to  a  single  surface.  When  multiple  surfaces  are  present 
in  each  pixel  or  detector,  the  process  of  accurately  identifying  the  correct  ranges  to 
each  surface  is  more  complex.  With  some  modification,  traditional  ranging  techniques 
can  account  for  multiple  surfaces  given  adequate  temporal  separation  of  the  received 
pulses.  However,  effects  due  to  the  medium  through  which  we  are  imaging  introduce 
error  into  a  traditional  ranging  technique.  By  diffracting  light  from  neighboring  areas 
within  the  scene,  each  detector  or  pixel  may  receive  a  pulse  in  which  numerous  false 
surfaces  can  be  identified  in  addition  to  the  potential  for  multiple  true  surfaces. 

This  chapter  is  organized  as  follows:  Section  4.1  details  the  techniques  employed 
to  find  solutions  for  the  components  of  the  Gaussian  mixture  model,  in  Section  4.2 
an  alternative  solution  is  derived  with  a  constraint  applied  to  the  amplitudes  of  the 
estimated  pulses,  Section  4.3  details  the  parameters  used  for  simulation,  Section  4.4 
presents  considerations  and  challenges  associated  with  employing  EM  techniques  to 
solve  this  problem  and  Section  4.5  presents  a  comparison  of  range  RMSE  for  various 
competing  algorithms  or  techniques. 
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4.1  EM  Solution 


The  following  algorithm  is  derived  using  an  EM  approach  to  jointly  estimate  the 
pulse  parameters  for  a  total  of  N  surface  returns  as  well  as  a  signal  bias.  The  MSID 
approach  is  similar  to  the  EM  approach  derived  by  Dolce  for  fusing  2-D  and  3-D 
LADAR  data  [16].  However,  the  work  is  not  equivalent  because  Dolce’s  work  does 
not  account  for  the  possibility  of  multiple  surfaces  per  pixel,  nor  does  it  provide 
an  estimate  for  the  pulse  amplitude  or  width.  The  EM  process  is  generally  broken 
down  into  two  steps.  First,  the  E-Step  involves  formulating  a  statistical  relationship 
between  the  data  collected  (incomplete-data)  and  some  known  data  model  (complete- 
data),  and  then  finding  the  expectation  of  the  complete-data  log- likelihood.  The 
second  step,  the  M-Step  is  to  iteratively  maximize  the  expectation  of  the  complete- 
data  log- likelihood  function  found  in  the  E-Step  [14].  The  following  sections  will 
describe  in  detail  the  process  for  jointly  estimating  the  pulse  amplitudes,  widths, 
ranges  and  overall  bias  for  each  image  pixel. 

4.1.1  Formulating  the  Complete  and  Incomplete-data. 

As  previously  mentioned,  the  first  step  to  the  EM  algorithm  is  to  formulate  the 
statistical  relationship  between  the  observed  or  collected  data  (incomplete-data)  and 
the  data  model  (complete-data).  Since  we  are  assuming  our  received  pulse  is  essen¬ 
tially  a  Gaussian  mixture  with  an  additive  bias,  the  relationship  is 

N  M  M 

d  ( u ,  v,  rk)  =  ^2  ^2  ^2  ^ n)  (“’ v ’  X1  y-> +  ^ B  (“’  v> rfc)  (4.1) 

n=  1  x=l  y=  1 

where  d  is  the  incomplete-data  and  the  complete-data  is  made  up  of  a  component 
for  the  received  pulse,  d^n\  and  a  separate  component  for  the  detector  bias,  ds-  It  is 
important  to  note  at  this  point  that  we  are  finding  a  total  of  N  amplitudes,  ranges  and 
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corresponding  pulse  widths  for  each  pixel.  The  expected  values  of  the  complete-data 
components  are  shown  in  (4.2)  and  (4.3). 
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d(n)  (u,  v,  x ,  y,  rk) 


4(n)  ( x ,  y)  - 
(x,  y) 


rfc— ( x,y ) 
2cr(n)  (x,y)^ 


x,  v  -  y) 


(4.2) 
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dB  (u,v,  rk) 


B  ( u ,  v ) 


(4.3) 


Since  both  components  of  the  complete-data  are  approximated  by  the  Poisson  PMF, 
and  the  sum  of  Poisson  random  variables  is  also  Poisson,  we  can  now  state  the  joint 
probability  of  the  complete-data,  pj,  as 


Vj 

x 


N  M  M  K  M  M 

=  |  n  n  n  n  n  n? 

n=l  u=  1  v=l  k= 1  x=l  y=  1 
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M  M  I\ 

n  n  U  pb 


\U=1  V=1 


k= 1 


\u,v,x,y, 


where  the  probability  for  the  nth  pulse,  p^n\  is 


(4.4) 


rk—  r^n)  ( x,y ) 

(  _  A^)  (x  v)  e  2aM(x,y)2  h(u-X,V-yf(n)  (u’v’*’y’rk) 

{din)  (W>  V’  r '*)  }  =  - i(»Hu,v,x,y,rk)l - 


XexP  {  y/£v&(x,y)e  2CT(n,(^)2  h{u-x,v-y)\  , 


rk—  r(n)  (a;, y) 


(4.5) 


and  the  probability  for  the  bias,  pB,  is 


Pb 


dB  (■ u ,  u,  r fe) 


dB  {u,v,rk)\ 


(4.6) 
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Now  that  we  have  the  joint  probability,  we  can  form  the  complete-data  log-likelihood, 
L,  by  taking  the  natural  logarithm  of  equation  (4.4)  such  that 


N  M  M  I<  M  M  , 

(u,v,x,y,rk) 

n=  1  u=  1  v=l  k=  1  x=l  y=  1  ' 


|  rk-r(n\x,y) 


x  ln  \  2a(n>(x’y)2  h  (u-X,V~y) 


,.-r(nh 


(4.7) 


Vsl’l)(h)e  h  (u  -  x,  v-y)- in \d^(u,v,x,  y,  rt)\ 


M  M  K 

+  EEE  \  dB  (u,  V,  rk)  ln  [B  (u,  u)]  -B(u,v)  -  In  |  dB(u,v,rk)\ 

u=  1  /c=l 


With  the  complete-data,  log-likelihood  formed,  we  are  now  ready  to  perform  the  E- 
Step. 


4.1.2  Finding  the  Expectation  E-Step. 

The  E-Step  involves  finding  the  expectation  of  (4.7)  conditioned  on  the  incomplete- 
data  and  the  previous  estimates  for  the  pulse  and  bias.  Through  the  course  of  the 
following  derivation,  it  was  realized  that  maximizing  the  expectation  with  respect 
to  the  bias  and  the  pulse  are  separable.  Therefore,  the  following  equations  will  be 
broken  up  into  respective  conditional  expectations  for  the  Gaussian  pulse  mixture, 
Q(n\  and  the  bias,  Qb,  to  simplify  the  explanation.  We  now  let  Q  be  the  overall 
conditional  expectation  of  the  complete-data  log- likelihood  function  such  that 

N 

q  =  Q(n) + «*-  <4-8) 

n=  1 
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where 


M  M  K  M  M 


Q(n)  =  £  £  £  E  E  [E  d(n)  (u,v,x,y,rk)  d  (u,  v,  rk) ,  P$}  (u,  v,  rk) ,  Bold  («,  v) 
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x  <  In 
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-E  I  In 


d(n)  (u,  v,x,y,rt)\  d  (u,  v,  rk) ,  (u,v,rk) ,  BoU  (u,v 


(4.9) 


and 


M  M  K  r 

Qb  =  E  E  E  E  dB(u,  v,  rk)  d  (u,  v,  rk) ,  P^J  («,  v,  rk) ,  Bold  (u,  v )  |  In  [5  («,  v)] 

n=l  i;=l  /c=l  L 


B  («,  v)  -  i?  ■!  In  dB(u,v,rk)\  d(u,v,rk) ,  P^J  (u,v,rk) ,  Bkld(u, 


(4.10) 


It  can  be  shown  that  the  expected  value  of  the  complete-data  components  with  re¬ 
spect  to  the  incomplete-data  and  the  previous  estimates  for  the  pulse  and  bias  are  as 
indicated  in  (4.11)  and  (4.12)  [62], 
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d(n)  («,  v,  x,  y,  rk)  |d  (u,  v,  rk) ,  (u,  v,  rk) ,  Boid  (■ u ,  v) 

d(u,v,rk)P{0u(x,y,rk)h(u-x,v-y) 
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(4.11) 
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dB  {u,  v,  rk)  | d  (u,  v,  rk) ,  (u,  v,  rk) ,  Bold  (u,  v 


d  (u,  v,r fc)  Bold  (u,v) 
I 0id  (■ u,v,rk ) 


(4.12) 


In  (4.11)  and  (4.12),  I0id(u,v,rk)  is  the  image  produced  by  the  pulse  estimate  and 
the  additive  bias  and  is  equal  to 


N  M  M 

hid  (u,  v ,  rk)  =  EEE  pdd  (x>  y,  rk)  h  (u  —  x,v  —  y)  +  Bold  (u,  v ) .  (4. 13) 

n=  1  x=l  y=  1 
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4.1.3  Maximizing  the  Expectation  M-Step. 


Now  that  we  have  the  expectation  of  the  complete-data  log- likelihood  function, 
we  can  maximize  it  with  respect  to  N  total  ranges,  pulse  widths  and  amplitudes  as 
well  as  the  signal  bias.  Similar  to  the  previous  work  by  Dolce,  the  received  pulse 
from  each  surface  is  assumed  to  exist  entirely  in  the  range  gate  [16].  This  assumption 
allows  us  to  find  a  direct  solution  for  range.  In  order  to  find  estimates  for  A^n\  An\ 
cr-D  and  B  we  must  take  the  derivative  of  Q  with  respect  to  each  parameter,  set  the 
resultant  derivative  equal  to  zero  and  then  solve  for  the  desired  parameter. 

Even  though  we  are  looking  for  a  solution  that  involves  N  surface  returns,  it 
is  possible  and  perhaps  even  likely  depending  on  the  imaged  scene  that  the  FOV 
for  a  given  pixel  may  have  fewer  surfaces  visible.  In  those  cases,  the  corresponding 
amplitudes  for  ranges  that  are  not  truly  present  in  the  image  scene  are  driven  towards 
zero  by  this  algorithm  as  will  be  demonstrated  in  Section  4.5.  Several  techniques 
exist  for  finding  an  upper  bound  on  the  number  of  surfaces  in  each  pixel  such  as  the 
center-of-gravity  and  zero-crossing  of  the  first  derivative  of  the  received  pulse  [67]. 
This  work  assumed  a  fixed  cap  of  two  visible  surfaces  per  pixel;  however,  this  cap 
could  be  adjusted  to  account  for  varying  numbers  of  surfaces  per  pixel  depending  on 
the  number  of  2-D  slices  per  image. 

We  will  first  demonstrate  the  maximization  process  for  the  pulse  amplitudes. 
Bringing  the  derivative  inside  of  the  summations  of  (4.14)  is  the  first  step  to  maximiz¬ 
ing  the  expectation  with  respect  to  amplitude.  Upon  inspection  it  is  evident  that  Qb 
in  (4.8)  is  not  dependent  on  (. x,y ),  thus  its  derivative  with  respect  to  A^  (. x,y ) 
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will  be  zero  and  we  are  left  with 
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(4.14) 


Computing  the  derivative  of  (4.14)  and  eliminating  terms  with  no  dependency  on 
(. x,y ),  we  can  further  simplify  to 


_  N  M  M  K  M  M  /  ,,  ,„(«),  \ur  \ 

8Q  =  ££££££(  d(W,U’rk'P°ld  ( x,y,rk)h(u-x,v-y )  ^ 


dA (n)  ( x,y ) 


n=l  w=l  u=l  /c=l  x=l  y=l 
\rk~r^  (x,y)  1 


Iold{u,v,rk) 


1  dA(n\x,y) 
A(n^(x,y)  dA(n^ (xo,yo) 


.__L_e  h(u-x,v-y) 


dA(n\x,y) 
dA(n)  (x0,yo) 


(4.15) 


We  now  have  that  the  derivative  of  Q  with  respect  to  ( x ,  y )  is 
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2CTW(x.y,2  h(u~x,v^y)  V  <J(x-x0,j/^j/0). 


The  Dirac  delta  function  in  (4.16)  allows  us  to  remove  two  of  the  summations  with 
respect  to  x  and  y  via  the  sifting  property,  since  all  values  not  equal  to  xq  and  ijq  will 
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be  zero,  leaving  us  with 
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We  can  now  set  (4.17)  equal  to  zero  and  solve  for  A ^  (, x,y ), 
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If  we  now  recall  the  assumption  that  the  pulse  is  entirely  within  the  range  gate  and 
make  the  additional  assumption  that 
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^2  h(x,y)  =  l,  (4.19) 

x,y=  l 


the  denominator  of  (4.18)  will  sum  to  one  leaving  us  with  a  final  iterative  solution 
for  A (x,  y)  that  is 
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We  can  perform  a  very  similar  process  for  maximizing  Q  with  respect  to  range,  but 
this  time  we  consider  our  assumption  of  the  pulse  being  within  the  range  gate  up 
front.  This  assumption  forces  the  summation  of  the  pulse  with  respect  to  a  change 
in  range  to  remain  a  constant.  Or  in  other  words,  even  though  the  pulse  location  is 
dependent  on  range,  the  summation  of  the  derivative  with  respect  to  r^n\xo,yo)  is 


zero.  We  now  just  have  the  derivative 
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to  maximize.  Following  a  very  similar  simplification  process  to  the  above  process  for 
the  amplitude,  we  are  left  with  the  solution  for  range 
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This  solution  is  nearly  identical  to  the  solution  found  by  Dolce,  with  the  exception 
that  it  is  generalized  for  multiple  surfaces  per  detector,  and  it  does  not  consider  the 
effects  of  3-D  and  2-D  fusion  [16].  Originally,  pulse  width  was  not  of  major  concern 
because  the  simulation  results  were  not  sensitive  to  changes  in  pulse  width  given 
that  all  targets  were  simulated  to  be  oriented  normally  to  the  illumination  source. 
However,  through  the  course  of  employing  this  algorithm  on  experimentally  collected 
data  where  targets  were  not  always  oriented  normally  to  the  illumination  source, 
it  was  realized  that  pulse  width  deviation  caused  a  noticeable  error  in  the  results. 
For  this  reason,  the  solution  for  pulse  width  was  derived  using  the  same  techniques 
previously  presented,  and  found  to  equal 
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(4.23) 


During  the  iterative  process,  the  algorithm  will  then  allow  the  pulse  width  to  adjust 
based  on  the  orientation  and  physical  properties  of  the  targeted  surfaces. 
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We  can  now  perform  a  similar  procedure  to  find  the  solution  for  pixel  bias.  For 
this  solution  we  again  refer  back  to  equations  (4.8),  (4.9)  and  (4.10).  We  notice  that 
(4.9)  is  not  dependent  on  pixel  bias.  Therefore  taking  the  derivative  of  (4.10)  with 
respect  to  B(u0,v0 )  will  be  sufficient 
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Equation  (4.24)  simplifies  to 
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The  partial  derivative  of  B(u,v )  with  respect  to  B(uq,v o)  will  yield  a  Dirac  delta 


allowing  for  further  simplification  to 
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This  equation  can  now  be  set  equal  to  zero  and  solved  for  B(uo,v o)  yielding  the 
solution  for  pixel  bias  as  shown  in  (4.27) 
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Using  the  results  from  (4.20),  (4.22),  (4.23)  and  (4.27)  we  are  now  able  to  iteratively 
ford  estimates  for  all  parameters  simultaneously.  The  following  section  presents  an 
added  constraint  on  amplitude  estimation  that  further  enhances  the  capability  of  the 
algorithm. 

Initialization  is  a  commonly  cited  challenge  for  EM  algorithms  [69].  However,  the 
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iterative  solutions  presented  in  (4.20),  (4.22),  (4.23)  and  (4.27),  can  effectively  be 
initialized  as  follows: 

•  Initialize  range  through  peak  search  or  correlation  technique 

•  Initialize  amplitude  to  intensity  values  at  initial  range  estimates 

•  Initialize  pulse  width  to  the  outgoing  pulse  width 

•  Initialize  pixel  bias  to  minimum  value  in  received  pulse 

This  methodology  for  initialization  proved  effective  as  shown  in  the  later  results. 

4.2  Constrained  Amplitude  EM  Solution 

It  is  generally  accepted  that  the  Richardson-Lucy  algorithm  (3.2)  will  produce  a 
maximum  likelihood  estimate  of  a  2-D  scene  when  the  blurring  source  or  PSF  is  known 
and  the  noise  is  Poisson  [57],  [62],  However,  when  an  image  is  both  spatially  and 
temporally  blurred  as  is  the  case  with  3-D  imaging  through  a  turbulent  atmosphere, 
the  Richardson-Lucy  algorithm  may  be  less  than  optimal. 

The  solutions  provided  in  Section  4.1.3  allowed  us  to  directly  solve  for  the  am¬ 
plitude,  range,  pulse  width  and  bias  for  N  surfaces  as  detected  by  each  pixel  in  a 
APD  array  while  also  removing  the  spatial  and  temporal  blurring  associated  with  the 
PSF.  Using  the  stopping  criteria  from  Section  2.4,  it  was  observed  that  the  algorithm 
often  converged  before  optimal  estimates  on  pulse  amplitude  were  achieved.  While 
the  results  from  the  algorithm  were  still  an  improvement  over  traditional  techniques, 
an  attempt  was  made  to  improve  upon  this  solution  by  considering  the  following 
constraint.  When  dealing  with  full- waveform  data,  each  3-D  image  can  be  flattened 
into  an  amplitude  only  2-D  image  by  simply  summing  the  amplitude  for  each  slice 
in  the  data  cube.  With  this  in  mind,  the  amplitude  only  2-D  image  obtained  from 
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the  Richardson- Lucy  deconvolution  algorithm  appears  to  be  a  potential  candidate 
for  an  added  constraint  to  prevent  inappropriate  convergence  of  the  EM  Algorithm 
designed  for  3-D  data.  Alternatively,  a  more  time  consuming  but  likely  more  accu¬ 
rate  approach  would  be  to  use  the  ML  estimate  derived  by  MacDonald  on  each  2-D 
frame  as  described  in  Section  3.4  prior  to  summation  into  an  overall  2-D  image  of  the 
remote  scene  [44],  Finally,  a  constraint  of  this  nature  could  also  be  applied  in  cases 
where  diffraction  is  not  of  concern.  In  this  case,  the  2-D  image  used  as  the  constraint 
would  be  obtained  by  summing  the  received  data  along  the  temporal  axis.  The  work 
by  Schulz  supports  the  idea  of  introducing  a  penalty  function  via  a  Lagrange  multi¬ 
plier  to  prevent  the  algorithm’s  convergence  to  undesirable  solutions  [60].  Provided 
the  function  that  we  add  as  a  penalty  is  continuously  differentiable,  the  method  of 
using  a  Lagrange  multiplier  will  allow  for  convergence  to  a  maximal  solution  for  the 
log-likelihood  subject  to  the  constraint  [65]. 

Using  the  amplitude  only  representation  of  the  3-D  data  cube,  we  found  the  best 
results  when  using  the  ML  estimate  on  each  2-D  frame  prior  to  summation  into  the 
constraint  image.  This  constraint  image  is  used  to  penalize  estimates  for  amplitude 
any  time  the  following  equation  is  not  satisfied 

N 

Ac  (x,  y)  =  A(n)  (x,  y ),  (4.28) 

n= 1 

where  Ac  is  the  amplitude  of  the  2-D  constraint  image.  The  original  multi-surface 
EM  algorithm  allowed  for  maximization  of  the  log-likelihood,  L ,  with  respect  to  N 
amplitudes,  ranges  and  pulse  widths  as  well  as  a  pixel  bias.  In  the  case  of  the  N 
surface  model,  this  can  be  written  in  simplified  notation  as 

N 

L  (A,  r,  a,B)=J2f  A”*.  r(”>,  B) ,  (4.29) 

71=1 
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where  /  represents  the  component  of  the  log-likelihood  for  each  of  the  surface  returns 
and  is  a  function  of  pulse  amplitude,  range,  width  and  bias.  We  now  want  to  maximize 
L  (A,  r,  cr,  B )  subject  to  a  new  constraint,  g  (A),  where 

N 

g(A)  =  J2AH-Ac  =  0.  (4.30) 

n=  1 

In  order  to  incorporate  this  constraint,  we  must  introduce  a  new  variable,  x,  f°r  our 
Lagrange  multiplier  and  then  maximize  the  constrained  log-likelihood,  L^, 

Lv  ( A ,  r,  a,  B,x)  =  L  ( A ,  r,  a,  B)  —  x  [g  (A)] .  (4.31) 


By  incorporating  the  constraint  via  the  method  above,  we  ensure  that  mathemati¬ 
cal  rigor  is  retained,  and  can  demonstrate  the  enhanced  performance  over  the  non- 
constrained  algorithm.  Using  (4.7)  as  a  basis,  Lv  becomes 
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With  the  constrained  complete-data  log-likelihood  formed,  we  are  again  ready  to 
perform  the  EM  process. 

It  is  readily  apparent  from  (4.32)  that  the  solutions  for  range,  pulse  width  and 
bias  will  not  be  changed  by  the  inclusion  of  the  constraint.  Therefore,  only  a  new 
solution  for  amplitude  and  the  constraint  parameter  are  needed.  Due  to  the  addition 
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of  the  constraint,  Q^n\  from  (4.9)  becomes 
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Ultimately  through  the  maximization  process,  we  find  that  the  constrained  solution 
for  amplitude  becomes 
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(4.34) 


If  we  now  recall  the  assumption  that  the  pulse  is  entirely  within  the  range  gate  as 


well  as  the  assumption  on  the  PSF  as  shown  in  (4.19),  the  right  hand  term  in  the 
denominator  of  (4.34)  will  sum  to  one  leaving  us  with  a  final  iterative  solution 
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(4.35) 


At  this  point,  the  only  thing  remaining  is  to  derive  the  solution  for  the  Lagrange 
multiplier.  The  derivative  of  Q ^  with  respect  to  x(xo>I/o)  easily  simplifies  to 
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(4.36) 
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Setting  (4.36)  to  zero  yields  the  following  relationship: 


N 

^2  Ain)  (®o,  Vo)  =  Ac  (xo,  Vo)  ■  (4.37) 

n=  1 


We  can  now  substitute  the  solution  from  (4.35)  into  this  relationship  to  allow  us  to 
find  a  solution  for  x(xo,Vo)- 
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The  resultant  solution  for  x(xo,Uo)  is 
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(4.39) 


At  this  point  we  have  all  of  the  necessary  pulse  parameter  solutions  in  both  a 
constrained  and  non-constrained  fashion.  By  incorporating  the  PSF  into  our  model 
we  have  enhanced  the  accuracy  of  the  model  and  enabled  the  removal  of  the  effects  of 
diffraction.  The  subsequent  results  presented  in  this  chapter  will  assume  that  the  PSF 
is  known.  However,  Chapter  V  will  demonstrate  that  the  average  PSF  parameterized 
by  ro  can  be  computed  with  the  help  of  this  algorithm. 


4.3  Simulation 

This  work  will  report  results  from  simulated  3-D  FLASH  LADAR  data  that  was 
designed  to  mimic  an  ASC  Tigereye™  camera.  For  this  reason,  simulation  parameters 
were  set  to  those  published  for  the  system. 
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4.3.1  Sensor  Parameters. 


Table  4.1  identifies  the  known  specifications  for  the  ASC  sensor.  Proper  sampling 

Table  4.1:  ASC  3-D  Tigereye  FLASH  LADAR  System  Specifications. 


Known  System  Parameters 

Parameter  Name 

Defined  Value 

Frames  per  image 

Laser  wavelength  (A) 

Sample  rate 

Range  delta  per  frame 

Total  range  gate 

Energy  per  pulse  ( Et ) 

Sensor  pulse  width 

Detector  Array  Size 

Pixel  Size 

20 

1.57  /jrn 

420  MHz 

0.357  m 

7.14  m 

0.025  J 
4.7x10-9s 

128  x  128 

100  /jrn  x  100  /jrn 

Lens  Parameters  for  Simulation 

Parameter  Name 

Defined  Value 

Sensor  focal  length  (/)) 
Aperture  diameter  ( D ) 
Instantaneous  Field  of  View  (iFOV) 

3  m 

2.325  cm 

0.24° 

of  the  images  is  important  due  to  the  use  of  deconvolution  in  the  algorithm.  The 
maximum  dimension,  7,  on  a  pixel  in  the  detector  array  must  abide  by  the  relationship 
[22] 

7  <  ^  (4.40) 

The  factor  of  two  in  the  denominator  arises  from  the  Nyquist  sampling  theorem 
which  states  we  must  sample  at  twice  the  maximum  achievable  spatial  frequency, 
ismax,  as  shown  in  (2.7).  Therefore,  given  the  known  system  specifications  for  the 
ASC  system,  the  lens  parameters  were  chosen  as  also  shown  in  Table  4.1.  While  this 
is  not  a  current  commercially  available  lens  configuration,  the  flexibility  of  simulation 
allows  for  selection  of  the  optics  parameters  in  order  to  ensure  proper  sampling  of 
the  simulated  images  according  to  (4.40).  Ultimately  the  goal  was  to  simulate  the 
type  of  diffraction  that  would  be  experienced  by  an  airborne  platform  incorporating 
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FLASH  LADAR  technology  as  a  remote  sensor  at  long  slant  ranges  from  the  target. 
Therefore  the  parameter  values  for  the  simulation  were  chosen  to  allow  the  ratio  of 
aperture  diameter,  D,  to  atmospheric  seeing  parameter,  ro,  to  be  similar  to  what 
would  be  experienced  in  an  airborne  sensor  application. 

4.3.2  Target  Profiles. 

All  targets  for  this  simulation  were  designed  such  that  all  surfaces  would  produce 
returns  within  the  range  gate. 
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(a)  3-bar  Target  Range  Image 


(b)  Obscured  Target  Range  Image 


(c)  Multiple  Void  Target  Range  Image 


Figure  4.1:  (a)  The  3-bar  target  has  an  opaque  background  at  a  distance  of  2  meters 

from  the  three  separate  opaque  raised  surfaces,  (b)  Has  a  foreground  at  a  distance  that 
is  2  meters  from  the  partially  obscured  surface  in  the  background.  Only  the  center 
section  of  the  foreground  is  transparent,  the  remainder  of  the  surfaces  are  opaque,  (c) 
Has  an  opaque  foreground  at  a  fixed  distance,  with  an  opaque  background  that  varies 
from  1  to  4  meters  from  the  foreground. 
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While  the  algorithm  has  demonstrated  increased  performance  in  multiple  surface 
ranging  for  all  tested  profiles,  this  dissertation  will  consider  three  unique  target  profiles 
shown  in  Figures  4.1  and  4.2.  The  target  profiles  were  designed  to  ensure  that  some 
pixels  had  various  numbers  of  surfaces  visible.  In  the  case  of  the  3-bar  target  and 
multi-void  target,  only  the  edges  of  the  features  possessed  multiple  surfaces  per  pixel. 
However,  the  obscured  target  has  an  entire  region  in  the  center  of  the  target  where 
multiple  surfaces  are  present.  For  the  example  targets  shown,  no  pixel  will  have  more 
than  two  surfaces  visible.  However,  once  the  effects  of  diffraction  are  added  in,  some 
pixels  may  receive  returns  from  additional  surfaces  depending  on  the  severity  of  the 
simulated  turbulence.  Given  the  relatively  short  range  gate  for  a  single  3-D  image, 
the  number  of  surfaces  in  the  included  target  profiles  was  limited  to  two.  However,  for 
a  larger  range  gate  the  algorithm  can  easily  be  expanded  to  account  for  the  general 
case  of  N  surfaces  [53]. 
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(a)  3-bar  Target  Range  Image 
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(c)  Obsured  Target  Range  Image 


X  Position 

(e)  Multiple  Void  Target  Range  Image 


X  Position 


(b)  Pixels  with  Multiple  Surfaces 


X  Position 


(d)  Pixels  with  Multiple  Surfaces 


X  Position 

(f)  Pixels  with  Multiple  Surfaces 


Figure  4.2:  (a)  Has  a  background  at  a  distance  of  2  meters  from  the  three  separate 

raised  surfaces,  (b)  The  areas  in  white  indicate  the  pixels  whose  iFOV  contain  2 
surfaces,  while  all  remaining  pixels  have  only  a  single  surface  visible  in  the  iFOV.  (c) 
Has  a  foreground  at  a  distance  that  is  2  meters  from  the  partially  obscured  surface  in 
the  background,  (d)  The  areas  in  white  indicate  the  pixels  where  the  iFOV  contains  2 
surfaces,  (e)  Has  a  foreground  at  a  fixed  distance,  with  a  background  that  varies  from 
1  to  4  meters  from  the  foreground,  (f)  The  areas  in  white  indicate  the  pixels  where 
the  iFOV  contains  2  surfaces. 
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4.4  Mixture  Modeling  Considerations 


Expectation  Maximization  in  conjunction  with  Gaussian  mixture  modeling  is  a 
common  technique  employed  to  isolate  the  parameters  of  interest  in  the  received 
LADAR  data  [46].  However,  Zhuang  points  out  that  common  challenges  with  employ¬ 
ing  EM  in  conjunction  with  Gaussian  mixture  modeling  are  determining  the  number 
of  components  in  the  mixture  or  isolating  mixture  components  that  may  merge  due  to 
their  individual  parameters  [69].  For  various  applications  of  Gaussian  mixture  mod¬ 
eling,  Vlassis  and  Nikas  solved  the  first  problem  with  a  Greedy  EM  approach.  They 
performed  an  iterative  process  where  the  number  of  components  was  incrementally 
increased  until  the  number  corresponding  with  the  solution  with  the  highest  likeli¬ 
hood  was  obtained  [66].  Unfortunately,  to  employ  this  solution  on  the  multi-surface 
detection  problem  including  the  effects  of  diffraction  would  be  intractable  due  to  the 
sheer  volume  of  possibilities.  For  instance,  if  we  consider  a  128  x  128  array,  where 
each  detector  may  have  between  0  —  2  surfaces  visible,  we  would  have  to  consider 
gi6384  poggibie  combinations  for  each  image  we  wish  to  process.  Additionally,  it  also 
seems  intuitive  that  the  maximum  likelihood  estimate  will  occur  where  the  maximum 
possible  number  of  surfaces  are  estimated  since  the  algorithm  will  attempt  to  fit  the 
noise  inherent  in  the  received  signal.  This  research  proposes  a  new  solution  to  this 
problem.  The  Cap  and  Refine  (CaR)  strategy  places  an  upper  bound  on  the  number 
of  surfaces  for  which  pulse  parameters  are  generated.  The  respective  pulse  amplitudes 
will  then  be  used  to  refine  the  number  of  surfaces  visible  in  each  detector. 

4.4.1  Key  Challenges  with  Multisurface  Modeling. 

In  practice,  the  true  upper  bound  on  the  number  of  surfaces  visible  by  each  detec¬ 
tor  could  be  established  through  various  methods.  Techniques  such  as  the  center-of- 
gravity  and  zero-crossing  of  the  first  derivative  have  been  employed  in  the  past  with 
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success  [67].  Based  on  the  objectives  of  this  research,  the  upper  cap  on  the  number  of 
surfaces  was  fixed  at  two,  but  could  easily  be  adjusted  to  account  for  the  possibility 
of  additional  visible  surfaces  per  pixel.  The  novelty  for  the  approach  employed  by 
this  research  does  not  arise  from  the  ability  to  define  the  upper  cap  on  the  number 
of  surfaces.  Rather,  the  novelty  is  that  we  can  use  the  unique  functionality  of  the 
MS1D  algorithm  in  conjunction  with  various  detection  schemes  to  efficiently  estimate 
the  number  of  components  in  our  Gaussian  mixture.  Or  in  other  words,  once  we  have 
solved  for  the  maximum  possible  components  or  surface  returns,  we  can  then  discard 
those  that  are  insignificant  after  a  single  execution  of  the  MSID  algorithm,  rather 
than  iteratively  search  for  the  best  possible  Gaussian  mixture. 


Range  (meters)  Range  (meters) 

(a)  Estimate  not  Accounting  for  Diffraction  (b)  Estimate  Accounting  for  Diffraction 

Figure  4.3:  Comparison  of  pulses  when  the  effects  of  diffraction  were  not  accounted  for 

(a) ,  and  when  the  effects  of  diffraction  were  incorporated  into  the  estimation  algorithm 

(b) . 

The  MSID  algorithm  produces  a  number  of  amplitude,  range  and  pulse  width 
estimates  for  each  detector  based  on  the  upper  cap  established.  In  cases  where  one  or 
more  of  the  surfaces  either  arises  from  effects  due  to  diffraction,  or  in  cases  where  there 
are  no  surfaces  visible  in  the  FOV,  the  corresponding  estimates  for  amplitude  will  be 
driven  towards  zero.  Due  to  the  effects  of  noise  and  the  residual  effects  from  diffraction 
even  after  deconvolution,  the  algorithm  may  not  perfectly  drive  undesired  surfaces  to 
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zero.  In  Figure  4.3  we  consider  two  scenarios  again  with  the  3-bar  target  from  Figure 
4.2(a).  In  Figure  4.3  (a)  and  (b),  we  plot  the  received  pulse  for  pixel  (55,75)  in  both 
an  ideal  case  where  the  effects  of  diffraction  are  not  present  (non-diffracted  pulse)  as 
well  as  the  actual  received  pulse  where  the  image  is  both  spatially  and  temporally 
blurred  due  to  the  effects  of  diffraction  with  -  w  2.  Since  the  selected  pixel  is  two 

r  o  1 


120- 

20  40  60  80  100  120 

X  Position 

Figure  4.4:  Identifies  the  location  of  the  image  pixel  used. 

pixels  above  the  edge  of  one  of  the  raised  bars  as  shown  in  Figure  4.4,  in  the  ideal 
scenario  we  would  only  expect  to  see  a  single  pulse  return  at  103  meters.  However,  the 
effects  of  diffraction  have  caused  some  of  the  light  that  reflected  off  of  the  background 
at  105  meters  to  fall  incident  onto  this  detector  as  well,  causing  a  temporal  distortion 
which  manifests  as  a  second  visible  return  in  the  received  waveform.  This  presents 
one  of  the  fundamental  problems  when  performing  multiple  surface  modeling  where 
the  collected  images  may  be  impacted  by  the  effects  of  diffraction. 

Figure  4.3  also  demonstrates  the  utility  of  two  techniques  that  could  be  employed 
for  this  multi-surface  estimation  problem.  First  we  consider  a  traditional  Gaussian 
mixture  model,  Figure  4.3(a),  which  does  not  incorporate  the  effects  of  diffraction  in 
the  estimation  of  the  pulse.  In  this  case  we  see  that  the  technique  estimates  a  pulse 
that  closely  matches  the  received  pulse.  The  technique  also  extracts  the  amplitude 
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and  range  information  for  each  of  the  two  components  in  the  mixture.  However,  once 
we  incorporated  the  effects  of  diffraction  through  the  MSID  algorithm,  we  were  able 
to  drive  the  estimate  for  the  second  pulse  towards  zero  as  shown  in  Figure  4.3(b). 
While  the  second  return  is  barely  visible  at  this  point,  its  amplitude  is  not  perfectly 
zero,  highlighting  another  key  challenge  associated  with  the  multi-surface  estimation 
problem. 

At  this  point,  a  decision  is  required  as  to  whether  or  not  the  reflection  from  a  sec¬ 
ond  surface  represented  by  an  amplitude  of  353  photons  is  of  interest  or  not.  Clearly 
a  classical  detection  approach  could  be  employed  where  the  threshold  is  determined 
by  the  noise  variance.  For  instance,  Stilla  et  al.  compared  the  pulse  amplitude  to 
the  standard  deviation  of  the  background  noise.  In  cases  where  the  pulse  amplitude 
exceeded  three  times  the  standard  deviation  of  the  noise  for  at  least  the  duration  of 
the  transmitted  pulse,  their  initial  classification  was  that  a  pulse  was  present  [63]. 
They  then  used  waveform-stacking  to  refine  their  classification  of  whether  or  not  a 
surface  of  interest  was  truly  present.  While  this  detection  algorithm  is  efficient  to 
implement,  the  choice  of  threshold  is  somewhat  arbitrary  given  the  statistics  used  in 
our  model.  Additionally,  the  technique  does  not  take  advantage  of  the  fact  that  we 
are  using  full  waveform  data. 

4.4.2  Using  Probability  of  False  Alarm  as  Key  Metric. 

While  perhaps  not  as  efficient  to  implement  as  a  simple  threshold  detection  scheme, 
detection  methods  using  waveform  data  can  be  shown  to  be  extremely  accurate  [58] . 
Additionally,  the  performance  of  the  technique  can  be  adjusted  to  mission  specific 
roles  by  simply  adjusting  the  parameter  of  interest.  For  instance  if  we  primarily 
wanted  to  minimize  the  possibility  of  detecting  a  surface  that  doesn’t  truly  exist,  we 
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could  adjust  the  probability  of  false  alarm,  pF.  Given  the  following  hypotheses: 


H0  -  No  surface  is  present  in  this  pixel  at  this  range 
H i  -  Surface  is  present  in  this  pixel  at  this  range 

the  probability  of  false  alarm  is 


pF=p[H1\H0],  (4.41) 

At  this  point,  we  will  be  working  with  the  estimates  produced  by  the  MSID  algorithm 
in  conjunction  with  the  measured  data.  Ultimately,  the  goal  is  to  determine  the 
amplitude  threshold  that  will  ensure  we  exceed  the  user  defined  threshold  for  pF. 
Any  surfaces  with  amplitude  estimates  below  this  threshold  could  be  discarded  and 
classified  as  the  algorithm’s  attempt  to  fit  a  pulse  to  the  noise. 

In  order  to  execute  this  detection  strategy  we  also  need  to  establish  the  Likeli¬ 
hood  Ratio  Test  (LRT).  For  purposes  of  this  research  we  will  consider  equal  prior 
probabilities  and  equal  costs  allowing 

A  =  — — | — —  >  1  say  Hi.  Otherwise  say  H0  (4.42) 

P(D\H0) 

where  D  is  a  pulse  return  with  amplitude  corresponding  with  the  threshold  we  wish 
to  test  [58].  If  prior  knowledge  of  the  scene  is  available,  this  test  could  be  adjusted  to 
potentially  allow  for  better  performance.  Given  the  pulse  width  and  bias  estimates 
that  were  found  using  the  MSID  algorithm,  we  can  find  the  amplitude  threshold  that 
satisfies  the  user  defined  value  for  pF.  Unfortunately  due  to  the  fact  that  we  are  using 
full  waveform  data,  a  closed  form  solution  for  the  amplitude  threshold  does  not  exist. 
However,  Monte  Carlo  methods  could  be  employed  to  obtain  this  threshold  [58]. 
Using  the  width  parameter  generated  from  the  MSID  algorithm,  we  can  simulate 
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numerous  independent  noisy  waveforms.  Using  the  measured  or  estimated  bias  level, 
we  can  predict  the  noise  variance  in  the  signal.  The  amplitude  of  the  simulated  noisy 
return  can  be  gradually  increased.  For  each  amplitude  level,  pp  is  computed  [58]. 
Our  amplitude  threshold  is  established  once  we  reach  the  user  defined  value  for  pF. 
In  this  manner,  we  can  compute  a  threshold  that  varies  with  the  width  of  the  received 
pulse. 

While  this  process  can  be  considered  computationally  intensive,  it  would  be  pos¬ 
sible  to  compute  the  threshold  in  advance  for  a  range  of  pulse  width  and  bias  values 
and  store  them  in  a  lookup  table  since  each  amplitude  threshold  is  entirely  dependent 
on  these  values.  Computing  the  values  in  advance  could  easily  allow  for  a  test  that 
could  be  executed  in  real  time.  Demonstration  of  this  technique  for  a  pp  =  0.01  in 
conjunction  with  the  performance  gains  from  the  MSID  algorithm  are  detailed  in  the 
following  section. 

4.5  Results 

The  following  results  will  demonstrate  the  performance  enhancements  of  the  algo¬ 
rithm  developed  through  this  research  over  traditional  Gaussian  mixture  modeling. 
Various  algorithms  were  tested  against  MCFA  images,  where  each  MCFA  image  is 
composed  of  30  separate  3-D  images.  Using  the  MCFA  images  serves  to  improve  the 
SNR.  First,  we  must  consider  a  metric  by  which  the  algorithms  will  be  compared. 

4.5.1  Range  Accuracy  Measurement. 

In  the  case  of  the  multi-surface  problem,  comparing  the  performance  of  the  tech¬ 
niques  is  not  simply  a  matter  of  looking  at  the  difference  between  the  range  estimate 
and  the  true  value.  We  must  also  consider  the  true  number  of  ranges  present  for  each 
pixel  and  the  algorithm’s  ability  to  accurately  predict  both  the  number  of  surfaces 
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visible  as  well  as  the  true  range  to  that  surface.  The  following  method  will  be  used  to 
compare  the  performance  of  the  algorithms.  For  this  work,  we  again  assume  that  the 
max  number  of  surfaces  visible  in  any  one  image  pixel  is  two.  However,  this  technique 
will  also  scale  as  the  number  of  surfaces  increases. 

In  a  3-D  FLASH  LADAR  system  we  are  primarily  concerned  with  the  accurate 
ranging  to  the  target,  but  amplitude  also  plays  a  key  role  in  visual  depiction  of  the 
target.  Since  the  algorithms  mentioned  above  are  initially  hard-coded  to  force  each 
pixel  to  have  the  maximum  number  of  surfaces,  we  must  balance  both  the  predicted 
range  and  amplitude.  The  following  range  accuracy  measurement  will  consider  the 
RMSE  between  the  predicted  and  true  ranges  weighted  by  the  predicted  amplitudes. 
In  that  manner  we  will  not  penalize  any  of  the  algorithms  for  predicting  a  false  range  if 
the  corresponding  amplitude  is  driven  to  zero  or  more  specifically  below  the  threshold 
established  by  the  techniques  listed  in  Section  4.4.2.  However,  if  an  algorithm  falsely 
predicts  a  range,  it  will  be  penalized  based  on  the  amplitude  corresponding  to  the 
falsely  predicted  range.  For  instance,  (4.43)  demonstrates  the  error  calculation  if  the 
algorithm  correctly  predicts  that  there  are  two  surfaces  present  for  a  pixel,  (4.44)  is 
the  error  calculation  for  a  pixel  that  only  has  one  surface  present  but  the  algorithm 
predicts  two  and  (4.46)  is  the  case  where  one  surface  is  present  and  the  algorithm 
correctly  predicts  this. 

e2,2  =  -A(1)  (x,y)  (r(1)  (x,y)  -  r{l){true)  (x,  y))2+A(2)  (x,y)  (r(2)  (x,y)  -  r(2)(true)  (. x,y ))2 

(4.43) 

62,1  =  -A(1)  (x,  y)  (r(1)  (x,  y)  -  r[true)  (. x ,  y))2  +  A(2)  (x,  y)  (r(2)  (x,  y)  -  r{true)  (. x ,  y))2 

(4.44) 

e1)2  =  A^  (x,y)  (r^  (x,y)  —  rR)hrae)  y))2+A(1^  (x,y)  (r(1)  (x,y)  —  r(2)(true)  ( x,y ))2 

(4.45) 

<q,i  =  (x,  y)  (r1^  (x,  y)  —  Atrue )  (x,  y))2  (4.46) 
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Once  the  error,  (predicted, true,  is  computed  for  each  pixel,  the  mean  error  is  taken  over 
the  entire  image  and  divided  by  the  mean  amplitude  of  all  estimated  surfaces.  Finally, 
the  square  root  of  this  value  is  taken,  giving  us  an  error  value  with  units  of  meters. 
This  final  value  will  be  used  to  judge  the  accuracy  of  the  algorithm. 

4.5.2  Comparison  of  Algorithms. 

This  section  will  compare  four  multi-surface  estimation  techniques  against  the 
three  target  profiles  listed  in  Figure  4.2.  The  primary  metric  that  will  be  used  to  judge 
overall  performance  will  be  the  range  accuracy  measurement  listed  in  Section  4.5.1. 
The  techniques  will  be  tested  against  images  with  various  levels  of  range  diversity  as 
well  as  two  levels  of  atmospheric  turbulence  based  on  ro  values  of  1  cm  and  2  cm. 


■  Gaussian  Mixture  ML  Frame-by-Frame 

M  Non-constrained _ M  Constrained _ 


3-bar  (1cm)  Obscured  (1  cm)  Multi-void  (1cm)  3-bar  (2cm)  Obscured  (2cm)  Multi-void  (2cm) 

Target  Type  with  Seeing  in  Parenthesis 


Figure  4.5:  Range  RMSE  Comparison. 

The  first  technique  under  consideration  was  a  simple  Gaussian  decomposition 
using  an  EM  technique.  This  technique  will  essentially  fit  the  best  possible  two- 
component  per  pixel  Gaussian  mixture  to  the  received  data  without  considering  the 
effects  of  diffraction.  The  second  technique  was  to  use  the  ML  estimate  (3.9)  on  each 
frame  of  the  3-D  image  prior  to  estimation  of  the  individual  surfaces.  Third,  we  will 
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look  at  the  performance  of  the  originally  developed  non-constrained  MSID  algorithm 
[53].  Finally,  we  will  look  at  the  constrained  MSID  algorithm  [52], 

From  the  results  in  Figure  4.5  it  is  visible  that  the  constrained  MSID  algorithm 
obtains  the  best  performance  for  the  estimation  problem  for  each  of  the  target  profiles 
considered.  The  overall  magnitude  for  RMSE  is  highly  dependent  on  target  type; 
however,  the  observed  trend  on  algorithm  performance  was  similar  across  all  target 
types.  The  improvements  of  the  constrained  MSID  algorithm  compared  with  non- 
constrained  MSID  algorithm  were  only  slight  in  some  cases.  However,  the  ability 
to  accurately  predict  the  number  of  surfaces  visible  by  the  detector  in  conjunction 
with  the  CaR  technique  was  significantly  enhanced  when  using  the  constrained  MSID 
algorithm. 

The  true  numbers  of  surfaces  visible  by  each  detector  under  ideal  conditions  are 
shown  in  Figure  4.2.  When  compared  with  the  truth,  the  constrained  MSID  algorithm 
clearly  outperforms  the  other  techniques  as  shown  in  Figure  4.6  for  the  3-bar  target 
blurred  by  an  OTF  with  an  ro  of  2  cm.  Here,  the  pixels  that  are  estimated  to 
have  two  surfaces  visible  are  indicated  in  white.  The  pixels  where  the  amplitude  of 
only  a  single  surface  is  determined  to  be  significant  are  either  gray  or  black.  The 
determining  factor  for  which  surface  is  used  to  estimate  an  individual  return  is  based 
on  the  initialization  value  for  range.  The  results  from  this  test  were  used  to  finalize 
the  optimal  mixture  of  components  based  on  the  estimated  data.  The  constrained 
MSID  algorithm  demonstrated  the  best  ability  to  drive  false  surfaces  towards  zero 
making  it  easier  to  develop  an  accurate  mixture  model  using  the  CaR  technique. 
For  comparison  purposes  similar  results  were  produced  for  the  obscured  target  in 
Figure  4.7  and  the  multi-void  target  in  Figure  4.8.  For  each  of  the  target  types,  the 
constrained  MSID  algorithm  demonstrated  the  best  ability  to  eliminate  false  surfaces 
that  appear  due  to  noise  or  diffraction  effects. 
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(a)  Gaussian  Mixture  (b)  ML  Frame-by- Frame 
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(c)  Non-constrained  MSID  (d)  Constrained  MSID 

Figure  4.6:  Surface  prediction  using  CaR  technique  with  r0  =  2 cm  for  (a)  Gaussian  mix¬ 
ture  EM  algorithm  without  including  the  effects  of  diffraction,  (b)  RL  Frame-by-Frame 
deconvolution  technique,  (c)  non-constrained  MSID  algorithm  and  (d)  the  constrained 
MSID  algorithm.  The  pixels  where  the  amplitudes  of  both  surfaces  yield  a  detection 
are  indicated  in  white.  The  pixels  where  the  amplitude  of  the  second  surface,  A^2>(x,y), 
was  the  only  one  of  significance  are  indicated  in  gray.  Finally,  the  pixels  where  the 
amplitude  of  the  first  surface,  AA>(x,y),  was  the  only  one  of  significance  are  indicated 
in  black. 
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(c)  Non-constrained  MSID 


(d)  Constrained  MSID 


Figure  4.7:  Surface  prediction  using  CaR  technique  with  r0  =  2 cm  for  (a)  Gaussian  mix¬ 
ture  EM  algorithm  without  including  the  effects  of  diffraction,  (b)  RL  Frame-by-Frame 
deconvolution  technique,  (c)  non-constrained  MSID  algorithm  and  (d)  the  constrained 
MSID  algorithm.  The  pixels  where  the  amplitudes  of  both  surfaces  yield  a  detection 
are  indicated  in  white.  The  pixels  where  the  amplitude  of  the  second  surface,  A^2>(x,y), 
was  the  only  one  of  significance  are  indicated  in  gray.  Finally,  the  pixels  where  the 
amplitude  of  the  first  surface,  AA>(x,y),  was  the  only  one  of  significance  are  indicated 
in  black. 
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X  Position 

(d)  Constrained  MSID 


Figure  4.8:  Surface  prediction  using  CaR  technique  with  r0  =  2 cm  for  (a)  Gaussian  mix¬ 
ture  EM  algorithm  without  including  the  effects  of  diffraction,  (b)  RL  Frame-by-Frame 
deconvolution  technique,  (c)  non-constrained  MSID  algorithm  and  (d)  the  constrained 
MSID  algorithm.  The  pixels  where  the  amplitudes  of  both  surfaces  yield  a  detection 
are  indicated  in  white.  The  pixels  where  the  amplitude  of  the  second  surface,  A^2>(x,y), 
was  the  only  one  of  significance  are  indicated  in  gray.  Finally,  the  pixels  where  the 
amplitude  of  the  first  surface,  AA>(x,y),  was  the  only  one  of  significance  are  indicated 
in  black. 
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Figure  4.9:  Comparison  of  3-D  surface  returns  for  3-bar  target  using  a  (a)  Gaussian 
mixture  EM  algorithm  without  including  the  effects  of  diffraction,  and  (b)  the  con¬ 
strained  MSID  algorithm  which  accounts  for  the  effects  of  diffraction.  By  including 
the  effects  of  diffraction,  our  estimate  of  the  3-D  target  is  significantly  more  accurate 
when  compared  to  the  truth  in  Figure  4.1  (a). 
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Figure  4.10:  Comparison  of  3-D  surface  returns  for  obscured  target  using  a  (a)  Gaus¬ 
sian  mixture  EM  algorithm  without  including  the  effects  of  diffraction,  and  (b)  the 
constrained  MSID  algorithm  which  accounts  for  the  effects  of  diffraction.  By  including 
the  effects  of  diffraction,  our  estimate  of  the  3-D  target  is  significantly  more  accurate 
when  compared  to  the  truth  in  Figure  4.1  (b). 


123 


(a)  Gaussian  Mixture 


(b)  Constrained  MSID 


Figure  4.11:  Comparison  of  3-D  surface  returns  for  multiple  void  target  using  a  (a) 
Gaussian  mixture  EM  algorithm  without  including  the  effects  of  diffraction,  and  (b)  the 
constrained  MSID  algorithm  which  accounts  for  the  effects  of  diffraction.  By  including 
the  effects  of  diffraction,  our  estimate  of  the  3-D  target  is  significantly  more  accurate 
when  compared  to  the  truth  in  Figure  4.1  (c). 
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In  Figures  4.12,  4.13  and  4.14  we  compare  the  range  error  by  pixel  for  the  various 
target  profiles  from  Figure  4.2  for  each  of  the  four  multi-surface  ranging  techniques. 
The  results  provide  an  intuitive  explanation  for  the  variation  in  range  RMSE.  The 
traditional  Gaussian  decomposition  strategy  has  minimal  range  error  when  there  are 
truly  two  surfaces  to  estimate,  such  as  in  the  center  of  the  image  for  the  obscured 
target.  However,  near  the  edges  of  the  obscuration  where  some  of  the  reflected  light 
from  the  far  surface  diffracts  onto  neighboring  pixels  we  experience  an  increase  in 
range  error  as  expected.  By  using  the  ML  estimate  on  each  frame  of  the  3-D  data, 
we  can  reduce  this  error  somewhat,  though  largest  reductions  in  error  occur  when  the 
effects  of  diffraction  are  directly  incorporated  into  the  multi-surface  ranging  algorithm. 
On  the  3-bar  target  and  the  multi- void  target,  we  again  see  that  as  the  light  from 
different  surfaces  is  diffracted  to  neighboring  areas  of  the  image,  we  experience  higher 
levels  of  error. 


Figure  4.12:  3-Bar  target  range  error  by  pixel  comparison  (meters2). 
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Figure  4.13:  Obscured  target  range  error  by  pixel  comparison  (meters2). 
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Figure  4.14: 


Multi- void  target  range  error  by  pixel  comparison  (meters2). 
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4.6  Chapter  Summary 


By  incorporating  the  effects  of  diffraction,  the  MSID  algorithm  developed  through 
this  research  allows  for  significant  enhancement  to  the  multi-surface  estimation  prob¬ 
lem  when  properly  sampled  images  are  taken  through  atmospheric  turbulence.  Two 
variations  of  the  MSID  algorithm  were  developed.  First,  direct  solutions  for  the  pulse 
return  parameters  and  pixel  bias  were  derived.  Second,  a  constraint  on  the  amplitude 
was  applied  which  allowed  for  much  more  refined  pulse  return  estimates.  Rather  than 
employ  deconvolution  techniques  that  are  tailored  for  2-D  images,  both  approaches 
incorporate  the  effects  of  diffraction  into  the  Gaussian  mixture  model.  Additionally, 
the  MSID  algorithm  simultaneously  solves  for  range,  pulse  width  and  amplitude  for 
multiple  surfaces  per  detector  while  enhancing  pulse  returns  that  may  have  been 
diminished  due  to  the  diffractive  effects  of  the  atmosphere. 

Through  the  incorporation  of  the  effects  of  diffraction  into  the  algorithm,  the  ill 
effects  of  temporal  and  spatial  distortion  were  simultaneously  reduced.  The  results 
obtained  from  the  MSID  algorithm  were  superior  to  those  obtained  through  more 
traditional  techniques.  Simulation  examples  show  that  the  MSID  algorithm  derived 
in  this  work  improves  range  estimation  over  standard  Gaussian  mixture  modeling 
and  frame-by-frame  deconvolution  on  average  by  89%  and  85%  respectively  based  on 
range  RMSE  calculations.  Current  limitations  on  the  technology  limit  the  ill  effects 
when  imaging  through  turbulence,  but  as  technology  improves  and  resolution  of  the 
detectors  increases,  the  ill  effects  will  become  significantly  more  pronounced. 

Given  the  technical  challenges  associated  with  the  manufacturing  of  these  sensors, 
the  issue  of  proper  sampling  appears  to  be  a  considerable  hurdle  to  overcome  for  many 
applications.  However,  the  goal  of  this  research  was  to  demonstrate  the  promise  of  the 
novel  technique  presented  in  this  chapter  on  multiple  surface  ranging  in  the  presence 
of  atmospheric  aberrations.  Additionally,  the  MSID  algorithm  would  have  current 
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applicability  where  sensors  with  very  long  focal  lengths  are  employed.  For  instance, 
employing  a  3-D  FLASH  LADAR  in  conjunction  with  an  astronomical  telescope  could 
facilitate  the  sampling  required.  An  optical  system  of  this  nature  could  yield  benefits 
in  SSA. 
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V.  3-D  FLASH  LADAR  Parameterized  Blind  Deconvolution 


This  chapter  focuses  on  the  development  of  a  novel  blind  deconvolution  algorithm 
employed  on  properly  sampled  3-D  FLASH  LADAR  images.  This  research  will  build 
upon  the  MSID  algorithm  previously  developed  to  minimize  the  effects  of  diffraction 
on  3-D  FLASH  LADAR  while  producing  accurate  ranging  to  multiple  surfaces  [52], 
Using  an  enhanced  version  of  this  algorithm  and  considering  the  range  diversity  in¬ 
herent  in  3-D  images  allows  for  simultaneous  estimation  of  the  parameterized  PSF 
and  spatial  /  temporal  enhancement  of  the  image.  Simulation  results  will  be  pre¬ 
sented  to  demonstrate  the  utility  of  this  algorithm  in  controlled  cases  where  the  true 
PSF  is  known  [51].  Experimental  results  will  also  be  presented  where  the  PSF  can 
be  measured. 

As  previously  mentioned,  parameterized  blind  deconvolution  for  2-D  images  is 
an  ill-posed  problem.  While  solutions  have  been  developed  that  can  overcome  this 
hurdle,  additional  assumptions  or  approximations  are  often  required.  Using  a  likeli¬ 
hood  maximization  approach,  this  research  will  show  that  by  adding  range  diversity 
through  3-D  imaging,  parameterized  blind  deconvolution  is  no  longer  an  ill-posed 
problem. 

This  chapter  is  organized  as  follows:  Section  5.1  details  the  system  of  equations 
that  leads  to  an  over-determined  problem  when  dealing  with  properly  sampled  FLASH 
3-D  images.  In  Section  5.2  the  strategy  for  finding  the  parameterized  OTF  is  provided. 
Section  5.3  presents  the  parameters  and  target  profiles  used  for  simulation  and  Section 
5.4  presents  the  findings  when  using  actual  experimentally  collected  images  impacted 
by  atmospheric  turbulence. 
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5.1  Joint  Estimation  of  Image  and  Atmospheric  Seeing  -  System  of  Equa¬ 
tions 

Previous  work  presented  in  Chapter  IV  assumed  that  the  PSF  was  known,  and 
simultaneously  provided  iterative  solutions  for  pulse  amplitude,  range,  width  and 
bias.  Using  the  MSID  algorithm  we  will  now  show  that  through  the  simultaneous 
estimation  of  all  of  the  visible  surfaces  for  each  detector,  we  can  also  accurately 
identify  the  PSF  parameterized  by  r0.  Both  the  constrained  and  non-constrained 
MSID  algorithm  have  been  tested  with  success.  However,  the  advantage  with  using 
the  non-constrained  algorithm  for  this  application  is  that  we  do  not  need  to  first 
perform  a  deconvolution  on  the  2-D  representation  of  the  received  image.  A  potential 
employment  strategy  would  be  to  first  use  the  non-constrained  MSID  algorithm  to  find 
the  optimal  parameterized  PSF.  We  could  then  use  this  optimal  parameterized  PSF 
in  conjunction  with  the  constrained  MSID  algorithm  to  further  refine  the  estimate. 

It  is  through  this  joint  estimation  that  we  are  able  to  fully  demonstrate  the  utility 
of  this  multiple  surface  ranging  capability  for  minimizing  the  spatial  and  temporal 
blurring  in  a  tactical  environment.  Using  the  average  atmospheric  models  in  (2.10) 
and  (2.11)  the  OTF  is  reduced  to  a  single  unknown,  r0.  The  ability  to  find  this  single 
unknown  for  the  OTF,  and  its  Fourier  relationship  to  the  PSF  can  be  considered  a 
solution  to  the  parameterized  blind  deconvolution  problem. 

While  the  problem  of  parameterized  blind  deconvolution  is  ill-posed  for  2-D  im¬ 
ages,  for  3-D  images  it  is  possible  that  the  problem  may  be  over-determined.  The 
over-determined  nature  of  the  problem  is  derived  from  the  fact  that  we  are  using 
3-D  images  where  the  total  image  is  collected  in  an  extremely  short  time  span.  The 
laser  pulse  is  of  such  a  short  duration,  that  the  atmosphere  can  be  considered  static 
and  a  single  PSF  is  applied  to  all  range  slices  in  the  image  [21],  The  2-D  system  of 
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equations  is 


M 

I(u,v)=  J2o(x,y)h(u-x,v-y\r0)  (5.1) 

x,y= 1 

where  i  is  the  image  estimate,  o  is  the  object  we  are  trying  to  estimate  and  our  PSF, 
h,  is  parameterized  by  r'o •  Here  we  have  a  total  of  M 2  equations,  but  we  have  (M2  +  l) 
unknowns,  since  for  every  image  we  have  both  an  object  sampled  by  M2  pixels  and 
a  unique  PSF.  The  3-D  system  of  equations  is 

N  M 

I  (u,  v,  r'k)  =  P ^  ( x ,  y,  rfc)  h  (u  —  x,v  —  y\  r0)  +  B  (u,  v )  for  k  —  {1  :  K} 

n= 1  x,y=l 

(5.2) 

where  K  is  the  total  number  of  range  samples  or  frames  in  the  image.  For  3-D  images 
we  therefore  have  KM 2  equations  but  only  ( 31V  +  1)M2  +  1  unknowns,  since  we  want 
to  estimate  amplitude,  range  and  pulse  width  for  each  return,  bias  for  the  detector 
and  ro  for  the  PSF.  Therefore,  we  now  have  the  possibility  for  an  over-determined 
problem  if  the  condition 

K  >  31V  +  2  (5.3) 

is  satisfied,  since  K  must  be  an  integer.  In  addition  to  the  condition  in  (5.3),  the 
targeted  scene  must  also  have  range  diversity  in  order  to  prevent  an  ill-posed  problem. 

5.2  Maximum  Likelihood  Solution  for  Atmospheric  Seeing 

The  goal  of  this  work  is  to  show  that  likelihood  can  be  maximized  through  joint 
estimation  of  range,  amplitude,  pulse  width,  bias  and  atmospheric  seeing.  Given  the 
joint  probability  of  the  received  data  as  shown  in  (2.6)  the  log-likelihood,  L ,  can  be 
computed  as 

M  K 

L=^2  (u,v,rk)hi(Itot(u,v,rk))  -  Itot  (u,v,rk)  -  In  (d  (u,  v,  rfc)!)].  (5.4) 

u,i 1=1  k= 1 


131 


In  (5.4),  the  final  term,  d(u,v,rk)\  is  a  constant  that  does  not  vary  as  we  search  for 
the  correct  value  of  tq  to  maximize  likelihood.  Therefore,  we  simply  seek  to  maximize 
this  adjusted  likelihood  function,  Lv, 

M  K 

[d  (u,  v,  rk)  In  (Itot  (u,  v,  rk))  -  Itot  (u,  v,  rk)}.  (5.5) 

u,v=l  k=  1 

5.3  Simulation 

This  work  will  report  results  from  simulated  3-D  FLASH  LADAR  data  that  was 
again  designed  to  mimic  an  ASC  Tigereye  camera.  For  this  reason,  simulation  pa¬ 
rameters  were  set  to  those  published  for  the  system  as  shown  in  Table  4.1.  This 
section  will  be  broken  down  into  three  main  parts.  First  we  will  present  the  target 
profiles  considered.  Second,  we  will  show  how  range  diversity  is  critical  to  the  ability 
to  jointly  estimate  the  pulse,  bias  and  atmospheric  seeing.  Finally,  we  will  show  that 
joint  estimation  of  multiple  surfaces  is  required  in  order  to  solve  this  problem. 

5.3.1  Simulated  Target  Profiles. 

Each  simulated  3-D  image  will  consist  of  20  individual  data  frames.  With  a 
sample  rate  of  420  MHz,  each  128  x  128  pixel  data  frame  will  represent  a  range 
delta  of  approximately  0.357  m  for  a  total  range  gate  of  7.14  m.  All  targets  for 
this  simulation  were  designed  such  that  the  surfaces  would  produce  returns  within 
the  range  gate.  Target  profiles  were  selected  to  illustrate  the  dependence  on  range 
diversity.  This  paper  will  consider  three  unique  target  profiles  shown  in  Fig.  5.1. 
With  the  exception  of  the  flat  target  in  Fig.  5.1(a),  the  target  profiles  were  designed 
to  ensure  that  there  was  range  diversity  throughout  the  scene.  The  flat  target  was 
designed  to  show  that  likelihood  cannot  be  maximized  using  the  techniques  described 
in  this  paper  without  range  diversity  in  the  scene. 


132 


.105 


(c)  Single  bar  target  -  reflectivity 


(d)  Single  bar  target  -  range 


(e)  Multi  void  target  -  reflectivity 


(f)  Multi  void  target  -  range 


Figure  5.1:  (a)  The  flat  target  has  a  single  bar  in  the  center  with  lower  reflectivity, 

(b)  The  flat  target  is  at  a  range  of  105  meters  across  the  entire  sensor  field  of  view,  (c) 
The  single  bar  target  again  has  a  single  bar  in  the  center  with  lower  reflectivity,  (d) 
Additionally,  the  single  bar  target  has  range  diversity  since  the  background  is  at  105 
meters,  and  the  single  raised  bar  in  the  center  is  at  103  meters,  (e)  The  multiple  void 
target  has  numerous  voids  with  lower  reflectivity,  (f)  The  multi  void  target  has  the 
entire  foreground  at  100  meters  and  various  size  voids  at  ranges  between  101  and  104 
meters. 
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5.3.2  Range  Diversity  and  Likelihood  Maximization. 


The  simulation  results  were  obtained  by  first  executing  the  multi-surface  ranging 
algorithm  on  an  image  ensemble  consisting  of  an  average  of  30  individual  3-D  images 
across  a  range  of  r0  values  from  0.001  m  to  0.025  m.  Image  averaging  was  used  to 
improve  SNR  and  the  frames  were  properly  registered  making  the  employment  of  the 
average  short  exposure  OTF  (2.10)  valid.  The  adjusted  likelihood  was  then  computed 
according  to  (5.5),  and  the  optimal  solution  was  chosen  as  the  one  that  maximized 
likelihood.  The  results  below  will  focus  on  the  ability  to  accurately  estimate  the 
value  for  r0,  since  the  results  in  Chapter  IV  demonstrated  the  capability  of  the  MSID 
algorithm  given  a  known  value  for  r0. 

Figure  5.2  shows  the  adjusted  likelihood  with  respect  to  r0  for  each  of  the  targets 
identified  in  Figure  5.1.  In  each  case,  the  simulated  image  for  the  target  was  developed 
using  1  cm  for  the  true  value  of  r0.  As  expected,  the  flat  target  did  not  allow  for 
maximization  of  r0  at  the  correct  value.  Rather,  likelihood  continued  to  increase  with 
increasing  r0.  On  the  other  hand,  both  the  single  bar  target  and  the  multi  void  target 
allowed  for  maximization  of  likelihood  at  the  correct  value  of  1  cm. 
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(b)  Single  bar  target  -  Likelihood  vs  r0 


x  1 0  ’* 


rQ  (meters) 

(c)  Single  bar  target  -  zoom  view 


(d)  Multi  void  target  -  Likelihood  vs  ro  (e)  Multi  void  target  -  zoom  view 


Figure  5.2:  (a)  The  flat  target  has  a  continuously  increasing  likelihood  with  respect  to 

increasing  levels  of  ro  (b)  The  single  bar  target  has  a  small  amount  of  range  variance 
throughout  the  image,  yet  likelihood  is  maximized,  (c)  Upon  zooming  in,  likelihood  is 
clearly  maximized  at  the  correct  value  of  1  cm  for  the  single  bar  target,  (d)  The  multi 
void  target  has  far  more  range  diversity,  and  likelihood  is  again  maximized,  (e)  Upon 
zooming  in,  likelihood  is  clearly  maximized  at  the  correct  value  of  1  cm  for  the  multi 
void  target. 
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5.3.3  Joint  Estimation  of  Multiple  Surfaces. 


It  should  be  noted  that  accurate  results  are  far  more  likely  for  this  joint  estimation 
problem  if  simultaneous  estimation  of  all  surfaces  is  accomplished.  Since  blurring 
occurs  both  temporally  and  spatially,  we  must  account  for  both  if  we  want  to  identify 
the  correct  blurring  function  using  this  algorithm.  The  single  bar  target  was  designed 
such  that  each  pixel  would  have  at  most  two  surfaces  in  its  FOV.  Additionally,  only 
the  pixels  with  an  FOV  that  contained  the  edges  of  the  bar  would  have  more  than 
a  single  surface  if  diffraction  effects  were  not  present.  By  constraining  the  algorithm 
so  that  we  solve  for  at  most  one  surface  per  pixel,  the  ability  to  predict  an  accurate 
level  for  r0  is  significantly  degraded.  This  is  likely  due  to  the  fact  that  the  system  of 
equations  in  (5.2)  no  longer  applies.  Figure  5.3  shows  that  the  maximum  likelihood 
value  does  not  occur  at  the  correct  value  of  r0  =  1  cm  when  we  only  estimate  a  single 
surface  per  pixel.  Instead,  likelihood  is  maximized  at  tq  =  0.5  cm. 
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Figure  5.3:  When  only  considering  a  single  surface  per  detector  for  the  single  bar 
target,  likelihood  maximization  produces  a  low  estimate  for  the  value  of  tq. 
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5.4  Experimental 


The  simulation  results  discussed  in  the  previous  section  demonstrated  high  levels 
of  precision  in  the  ability  to  estimate  the  multiple  surface  model  in  conjunction  with 
the  atmospheric  seeing  parameter.  While  the  simulation  was  designed  to  closely 
replicate  the  results  from  an  actual  3-D  FLASH  LADAR  sensor,  the  complexity  of 
the  problem  was  reduced  based  on  the  target  geometry.  For  the  simulations,  all 
targets  were  oriented  normal  to  the  sensor.  Therefore,  any  pulse  width  expansion 
could  be  considered  associated  with  the  effects  of  diffraction.  However,  experimental 
conditions  considered  targets  with  various  orientations.  In  this  case,  the  received  pulse 
width  could  vary  due  to  the  angle  at  which  the  outgoing  pulse  strikes  the  surface. 
The  experimental  results  demonstrate  that  the  algorithm  is  capable  of  separating  the 
cause  of  pulse  width  variation. 

5.4.1  Sensor  Parameters. 

The  sensor  used  for  this  research  was  a  modified  ASC  Portable  3-D  FLASH 
LADAR  Camera  Kit™.  The  sensor  used  all  of  the  standard  components  in  this 
commercially  available  sensor;  however,  they  were  oriented  into  a  different  configu¬ 
ration  for  the  LISAF  TPS  as  shown  in  Figure  2.4.  The  only  significant  differences 
between  the  specifications  listed  in  Table  4.1  used  for  simulation  and  the  sensor  used 
for  this  experiment  is  with  the  optics.  The  optics  parameters  for  the  experiment 
are  listed  in  Table  2.2.  Given  the  sampling  requirements  in  (4.40)  and  the  optical 
specifications  in  Table  2.2,  the  maximum  pixel  size  for  proper  sampling  should  be  2 
^111x2  fin l.  However,  the  actual  sensor  detector  had  a  much  greater  pixel  size  at  100 
jumxlOO  fan. 

Previous  efforts  with  similar  FLASH  LADAR  systems  addressed  the  sampling 
requirement  by  significantly  restricting  the  aperture  diameter  with  a  mask  [16],  [?]. 
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Figure  5.4:  Frequency  response  as  a  function  of  spatial  frequency  for  7’o  levels  between 
0.0005m  and  0.002m. 

For  the  size  of  the  individual  detectors  in  this  sensor,  the  aperture  would  have  needed 
to  be  approximately  2  mm  for  proper  sampling  with  this  technique.  By  restricting  the 
aperture  that  far  we  would  have  greatly  reduced  the  amount  of  light  gathered  thus 
negatively  impacting  the  SNR.  As  an  alternative,  we  chose  to  use  a  highly  turbulent 
atmosphere  to  reduce  the  sampling  requirement. 

The  atmosphere  can  essentially  be  treated  as  a  low  pass  filter.  The  more  turbulent 
the  atmosphere  becomes,  the  lower  the  sampling  requirement  will  become  as  well. 
Considering  the  relationships  in  (2.9)  and  (2.10)  we  can  show  how  the  maximally 
observed  spatial  frequency  is  decreased  as  turbulence  is  increased  or  in  other  words 
as  r o  is  decreased.  Given  the  pixel  pitch  of  the  ASC  sensor,  the  maximum  spatial 
frequency  satisfying  Nyquist  criteria  is  5,000  (1/m).  Yet  the  optical  specifications  in 
Table  2.2  in  conjunction  with  the  relationship  in  (2.7)  yield  a  sampling  requirement 
of  305.7  x  103  (1/m).  Therefore,  the  goal  was  to  find  an  r0  value  that  produced  a 
cutoff  frequency  below  5,000  (1/m).  As  shown  in  Figure  5.4,  the  frequency  cutoff  is 
reduced  in  conjunction  with  a  reduction  in  tq. 
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Based  on  the  available  sensor  parameters,  r0  levels  of  approximately  0.0015m  will 
yield  proper  sampling.  Atmospheric  seeing  levels  of  this  nature  would  be  difficult  to 
find  naturally.  Therefore,  an  extremely  turbulent  atmosphere  was  generated  in  front 
of  the  aperture  using  a  60,000  British  Thermal  Unit  (BTU)  heat  source.  This  heat 
source  consistently  created  ro  levels  between  0.001m  and  0.002m.  We  have  found 
that  the  MSID  algorithm  has  some  capability  to  deal  with  slightly  undersampled 
data.  However,  the  extent  to  which  the  data  can  be  undersampled  with  this  technique 
remains  to  be  proven  as  a  future  research  topic. 

Placing  this  source  of  turbulence  in  front  of  the  aperture  had  another  benefit  for 
this  experiment.  Since  this  sensor  was  originally  designed  to  be  placed  in  a  pod  on  an 
aircraft,  the  lens  was  focused  at  infinity.  With  the  ranges  used  for  this  experiment, 
the  lens  being  focused  at  infinity  resulted  in  a  blurring  of  the  images.  The  focus  error 
OTF  is  similar  to  the  atmospheric  OTF  at  lower  frequencies  [22] .  However,  there  are 
additional  high  frequency  components  in  the  focus  error  OTF.  By  imaging  through 
the  turbulent  atmosphere,  the  added  high  frequency  components  were  filtered  out. 
The  resultant  OTF  was  measured  to  be  very  close  to  the  theoretical  short  exposure 
OTF  as  shown  in  the  following  results. 

5.4.2  Experimental  Target  Profiles. 

The  images  collected  for  this  research  consist  of  two  separate  target  configura¬ 
tions.  Figure  5.5  displays  the  first  configuration  and  its  associated  range  profile. 
This  configuration  was  designed  such  that  ideally  only  the  edge  pixels  where  the  two 
sheets  overlap  would  have  multiple  surfaces.  The  second  target  configuration  shown 
in  Figure  5.6  is  similar  to  the  obscured  target  simulated  in  Chapter  IV.  With  this 
target  configuration  the  front  sheet  of  plywood  had  the  center  section  removed.  An 
aluminum  screen  was  then  placed  over  this  opening  to  allow  a  portion  of  the  light 
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to  pass  through  and  a  portion  to  be  reflected.  A  second  sheet  of  plywood  was  then 
placed  behind  the  first  sheet  to  provide  a  second  return. 

For  target  configuration  1,  the  two  sheets  of  plywood  were  separated  in  range  by 
approximately  2.6  m.  The  sheets  were  each  oriented  perpendicular  from  the  sensor. 
Additionally,  the  sheets  were  positioned  such  that  there  would  be  some  overlap  at 
the  center.  Given  the  3°  FOV  for  the  sensor  and  the  range  to  the  target,  each  pixel 
will  correspond  with  a  spatial  area  of  approximately  6.5  x  6.5  cm2.  Therefore,  we 
would  expect  only  the  column  of  pixels  in  the  detector  array  corresponding  with  the 
overlap  to  have  more  than  one  surface  in  its  FOV  in  an  ideal  environment.  However, 
the  following  results  show  that  this  is  not  the  case. 

For  target  configuration  2,  the  two  sheets  of  plywood  were  separated  in  range  by 
approximately  3.4m.  The  sheet  in  the  foreground  had  an  opening  in  the  center  that 
was  covered  with  an  aluminum  mesh  screen.  A  second  sheet  of  plywood  was  oriented 
behind  and  parallel  to  the  first  sheet.  This  configuration  allows  us  to  demonstrate 
the  ability  to  detect  an  object  through  an  obscuration.  Given  the  range  to  this  target 
and  the  3°  FOV  of  the  sensor,  each  image  pixel  will  correspond  with  an  area  of 
approximately  6x6  cm2. 
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(c)  3-D  View  of  Target  Configuration  1 


Figure  5.5:  (a)  Plywood  in  foreground  at  approximately  156.4m  overlaps  a  second  sheet 
of  plywood  placed  at  approximately  159m.  (b)  Range  profile  for  target  configuration 
1.  (c)  3-D  representation  of  target  configuration. 
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(b)  Range  to  Surface  1 


(c)  Range  to  Surface  2 
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(d)  3-D  View  of  Target  Configuration  2 


Figure  5.6:  (a)  Plywood  in  foreground  at  approximately  145.8m  contains  a  mesh  screen 
in  the  center  to  allow  ranging  to  the  second  sheet  of  plywood  placed  at  approximately 
149.2m.  (b)  Range  to  the  most  dominant  surface  in  each  pixel,  (c)  Range  to  second 
most  dominant  surface  in  each  pixel,  (d)  3-D  representation  of  target  configuration. 
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Frame  6  Frame  7  Frame  8  Frame  9 


Figure  5.7:  Individual  2-D  frames  for  target  configuration  #1. 
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Figure  5.8:  Individual  2-D  frames  for  target  configuration  #2. 
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Photons 


5.4.3  Atmospheric  Seeing  -  Truth  Measurement  Technique. 

Both  configurations  were  designed  such  that  the  true  level  of  r0  could  be  measured 
and  compared  with  the  estimates  obtained  using  the  MSID  algorithm.  The  true  seeing 
parameter  is  measured  using  the  technique  identified  in  Section  2.5.  To  implement 
this  measurement  technique,  we  can  consider  a  range  slice  that  only  has  one  sheet 
of  plywood  in  view.  From  this  2-D  slice,  we  can  calculate  the  step  response  by 
measuring  the  change  from  background  to  illuminated  surface.  For  example,  when 
only  considering  the  intensity  image  corresponding  with  a  range  of  159m  in  the  first 
target  configuration,  we  get  the  2-D  image  and  measured  frequency  response  as  shown 
in  Figure  5.9.  Candidate  OTFs  across  a  suitable  range  of  r0  values  can  then  be 
compared  against  this  measured  OTF  to  find  the  one  that  produces  the  MMSE. 
Notice  that  the  measured  frequency  response  contains  high  frequency  components 
that  are  absent  from  the  theoretical  candidate  OTFs.  This  is  likely  due  to  the  noise 
present  in  the  measured  step  response. 
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(c)  Measured  Frequency  Response 

Figure  5.9:  (a)  Step  target  obtained  by  only  considering  the  image  intensity  at  a 

range  of  159m.  (b)  Impulse  response  computed  by  first  finding  the  step  response  or 
vertical  change  in  intensity  and  then  taking  the  first  derivative  of  the  step  response, 
(c)  Measured  frequency  response  with  an  overlay  of  the  best  fit  short  exposure  OTF 
(r0=O. 0012m),  and  a  range  of  candidate  OTFs. 
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5.4.4  Experimental  Results. 


Using  the  previously  mentioned  target  configurations,  we  now  execute  the  MSID 
algorithm  for  a  range  of  ro  values.  The  search  window  for  ro  will  consist  of  levels 
from  0.0001m  to  0.01m  with  an  increment  of  0.0001m.  One  of  the  key  advantages 
with  this  algorithm  is  that  execution  of  the  MSID  algorithm  for  each  level  of  ro  is 
independent.  This  allows  for  easy  parallelization  or  efficient  search  patterns  such  as 
the  pattern  described  in  Figure  3.4.  Using  the  estimated  parameters,  we  can  compute 
the  adjusted  likelihood  according  to  (5.5).  As  an  example,  Figure  5.10  shows  the 
adjusted  likelihood  for  each  level  of  r0  within  this  search  window  for  Trial  #1. 


(a)  Adjusted  Likelihood  vs  ro  (b)  Zoom  View  of  Adjusted  Likelihood 

Figure  5.10:  (a)  Adjusted  likelihood  for  a  range  of  ?’o  values  from  0.0001m  to  0.01m. 

(b)  Maximization  of  likelihood  occurs  at  0.0012m. 

Figure  5.10  shows  that  maximization  of  the  adjusted  likelihood  occurs  at  ro  = 
0.0012m,  which  was  what  was  expected  based  on  the  measured  level  of  rQ.  Of  the 
four  experiments  conducted,  each  were  accurate  to  within  0.0002  meters  as  shown 
in  Table  5.1.  More  importantly,  the  results  of  the  experiments  were  consistent  with 
the  results  from  simulations  where  we  had  significantly  more  control  over  the  target 
geometry,  simulated  turbulence  and  optical  specifications.  The  spatial  and  temporal 
blurring  for  target  configuration  ^2  had  less  of  an  impact  than  what  was  observed 
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Table  5.1:  Comparison  of  Measured  and  Estimated  r0  Values. 


Target  Profile  #1  r0  Values  (meters) 

Trial 

Measured 

MSID  -  Estimated 

CoV  -  Estimated 

Frame  11 

Frame  19 

1 

0.0012 

0.0012 

0.0012 

0.0013 

2 

0.0013 

0.0014 

0.0013 

0.0013 

Targe 

t  Profile  7^2  r0  Values  (meters) 

Trial 

Measured 

MSID  -  Estimated 

CoV  -  Estimated 

Frame  12 

Frame  22 

3 

0.0019 

0.0018 

0.0017 

0.0018 

4 

0.0019 

0.0017 

0.0016 

0.0017 

for  target  configuration  #1  as  indicated  by  the  levels  of  r0  shown  in  Table  5.1.  The 
experiments  for  target  configuration  #1  were  conducted  on  6  February  2012  when  the 
ambient  air  temperature  was  below  40°F.  The  experiments  for  target  configuration  ^2 
were  conducted  on  3  October  2012  when  the  ambient  air  temperature  was  above  70°F. 
Given  the  lower  temperature  difference  between  the  ambient  air  and  the  output  of 
the  60,  000  BTU  heat  source,  we  expect  the  values  of  r0  to  be  higher  when  the  outside 
air  is  warmer  as  indicated.  The  next  section  will  visually  show  the  improvement  in 
estimation  capability  using  the  MSID  algorithm. 

The  results  from  the  CoV  method  detailed  in  Chapter  111  are  also  listed  in  Table 
5.1.  Based  on  the  results  from  simulation,  shown  in  Section  3.4,  we  hand  selected  two 
frames  from  each  3-D  image  to  perform  the  CoV  method.  Recalling  that  simulations 
show  improved  accuracy  in  higher  contrast  images,  we  chose  frames  11  and  19  for 
target  configuration  #1,  and  frames  12  and  22  for  target  configuration  #2.  The  results 
from  the  CoV  method  are  similar  to  the  estimates  obtained  through  maximization 
of  likelihood  in  conjunction  with  the  MSID  algorithm.  This  highlights  that  the  CoV 
technique  is  a  viable  backup  in  cases  where  abnormalities  in  pulse  shape  may  make 
maximization  of  likelihood  in  conjunction  with  the  MSID  algorithm  unreliable. 
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5.4.5  Observed  Performance  Enhancement. 


The  overall  goal  of  this  technique  is  to  simultaneously  find  the  pulse  parameters 
and  value  of  r$  that  maximizes  likelihood  as  demonstrated  in  Section  5.3.  For  com¬ 
parison  purposes  we  will  also  consider  an  EM  solution  that  does  not  account  for  the 
effects  of  diffraction  and  only  seeks  to  fit  the  best  possible  Gaussian  mixture  to  the 
received  data.  Range  RMSE  measurements  are  not  possible  since  we  do  not  have 
truth  range  information  to  each  point  in  the  scene.  However,  we  can  assess  the  algo¬ 
rithm’s  ability  to  accurately  predict  the  correct  number  of  surfaces  for  each  point  in 
the  scene  based  on  target  geometry  as  well  as  the  ability  to  accurately  predict  r0. 

The  MSID  algorithm  relies  on  the  CaR  technique  to  identify  the  true  number  of 
surfaces  per  pixel  [52],  The  target  configurations  and  the  sensor  range  gate  were  set 
up  such  that  a  maximum  of  two  returns  should  be  visible  in  each  detector.  Therefore, 
the  MSID  algorithm  will  estimate  a  Gaussian  mixture  with  two  returns  for  each 
pixel.  In  cases  where  a  pixel  has  fewer  returns  than  the  pre-defined  cap,  the  algorithm 
should  ideally  drive  those  amplitudes  towards  zero.  As  previously  demonstrated,  once 
the  algorithm’s  termination  criteria  is  achieved,  we  then  execute  the  full-waveform 
detection  strategy  to  determine  if  the  individual  pulse  amplitudes  warrant  a  detection 
[58].  Given  the  target  geometry  in  configuration  #1,  we  would  only  expect  the  single 
column  of  pixels  corresponding  with  the  overlap  of  the  two  sheets  of  plywood  to  have 
multiple  returns  when  the  effects  of  diffraction  are  removed.  For  target  configuration 
#2,  we  expect  the  center  region  of  the  plywood  and  the  edge  of  the  front  sheet  of 
plywood  to  have  multiple  returns. 

Figure  5.11  shows  that  the  overlap  region  between  the  two  sheets  of  plywood  is  as 
many  as  four  pixels  wide.  At  a  range  of  approximately  157  m,  this  would  mean  that 
the  overlap  region  is  up  to  26  cm  wide,  providing  a  significant  amount  of  ambiguity  to 
the  true  overlap  point.  In  Figure  5.12  each  of  the  received  pulses  and  pulse  estimates 
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Figure  5.11:  Target  Configuration  #1  -  Number  of  surfaces  detected  for  each  pixel 
without  accounting  for  the  effects  of  diffraction. 


are  compared  for  the  highlighted  region  from  Figure  5.11.  With  the  exception  of  pixels 
(60,69)  and  (65,69),  each  of  the  received  pulses  have  two  returns  corresponding  with 
the  surface  prediction.  However,  based  on  target  geometry,  we  would  only  expect  one 
of  the  received  pulses  to  contain  multiple  surface  returns  in  the  absence  of  diffraction 
effects. 
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(c)  Pixel  (65,69) 


Figure  5.12:  Individual  pulse  reconstructions  compared  with  received  pulses  for  the 
highlighted  region  from  Figure  5.11  without  accounting  for  the  effects  of  diffraction. 
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By  incorporating  the  effects  of  diffraction,  we  can  more  accurately  estimate  the 
number  of  surfaces  per  pixel.  Figure  5.13  shows  the  number  of  surfaces  that  are 
detected  per  pixel  when  the  value  of  r o  which  maximized  the  likelihood  (ro  =  0.0012m) 
is  used.  By  accounting  for  the  effects  of  diffraction,  the  region  of  overlap  is  reduced 
to  a  single  column  in  most  places.  In  Figure  5.14  each  of  the  received  pulses  and 
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Figure  5.13:  Target  Configuration  #1  -  Number  of  surfaces  detected  for  each  pixel 
when  accounting  for  the  effects  of  diffraction. 

pulse  estimates  are  compared  for  the  highlighted  region  from  Figure  5.13.  Through 
the  use  of  the  MSID  algorithm,  the  amplitude  for  the  second  surface  is  driven  below 
the  detection  threshold  for  all  but  pixel  (62,69). 
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Figure  5.14:  Individual  pulse  reconstructions  compared  with  received  pulses  for  the 
highlighted  region  from  Figure  5.13  when  accounting  for  the  effects  of  diffraction. 
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An  additional  metric  for  comparison  is  based  on  the  known  target  height.  The 
sheet  of  plywood  on  the  right  is  1.22  m  tall,  and  the  range  to  this  sheet  of  plywood 
is  approximately  157  m.  Given  this  geometry  and  the  sensor  parameters  we  would 
expect  the  sheet  to  only  require  a  vertical  range  of  19  -  20  pixels  in  the  detector  FOV. 
However,  when  the  straight  Gaussian  mixture  model  is  applied  to  the  received  data, 
the  target  is  measured  at  22  -  23  pixels  tall.  Once  the  MS1D  algorithm  is  applied, 
the  target  only  requires  a  vertical  range  of  19  -  21  pixels  in  the  detector  FOV.  This 
reduction  in  blurring  around  the  edges  of  the  target  can  be  observed  by  looking  at 
the  difference  between  Figures  5.11  and  5.13  as  shown  in  Figure  5.15. 


Figure  5.15:  Target  Configuration  #1  -  Difference  in  the  number  of  surfaces  detected 
for  each  pixel  when  not  accounting  for  the  effects  of  diffraction  minus  the  results  from 
the  MSID  algorithm. 


We  now  demonstrate  similar  results  for  target  configuration  #2.  When  a  Gaussian 
mixture  model  is  applied  to  the  received  data  without  accounting  for  the  effects  of 
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diffraction,  the  CaR  technique  reveals  multiple  surfaces  in  many  pixels  where  only 
a  single  surface  should  exist.  In  an  ideal  environment,  a  single  row  /  column  of 
pixels  should  have  multiple  surfaces  visible  corresponding  with  the  top  /  left  sides 
of  the  target  respectively.  Additionally,  the  center  section  of  the  target  should  also 
have  multiple  surfaces  visible.  However,  when  we  fail  to  account  for  the  effects  of 
diffraction,  we  again  see  that  the  spatial  and  temporal  blurring  around  the  edges  of 
the  target  results  in  numerous  false  detections  as  shown  in  Figure  5.16. 
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Figure  5.16:  Target  Configuration  #2  -  Number  of  surfaces  detected  for  each  pixel 
without  accounting  for  the  effects  of  diffraction. 

By  incorporating  the  effects  of  diffraction,  the  estimated  number  of  surfaces  per 
pixel  can  again  be  brought  more  in  line  with  the  expectation  based  on  target  geometry. 
Figure  5.17  shows  the  number  of  surfaces  that  are  detected  per  pixel  when  the  value 
of  r0  which  maximized  likelihood  (r0  =  0.0019m)  is  used.  By  accounting  for  the 
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effects  of  diffraction,  the  overlap  regions  on  the  left  and  top  of  the  target  are  reduced 
in  most  places.  Additionally,  there  is  a  reduction  in  blurring  around  the  edges  of  the 


Pixels 


Figure  5.17:  Target  Configuration  #2  -  Number  of  surfaces  detected  for  each  pixel 
when  accounting  for  the  effects  of  diffraction. 


target  as  observed  by  comparing  the  difference  in  Figures  5.16  and  5.17  as  shown  in 
Figure  5.18. 
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Figure  5.18:  Target  Configuration  #2  -  Difference  in  the  number  of  surfaces  detected 
for  each  pixel  when  not  accounting  for  the  effects  of  diffraction  minus  the  results  from 
the  MSID  algorithm. 
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5.4.6  Joint  Estimation  Requirement. 


The  previous  results  demonstrated  that  a  likelihood  maximization  approach  to 
the  parameterized  blind  deconvolution  problem  is  possible  with  properly  sampled  3- 
D  images.  As  a  final  consideration,  the  requirement  to  fully  account  for  all  surfaces 
visible  in  the  received  data  will  be  demonstrated.  Through  simulation,  we  found 
that  the  ability  to  accurately  estimate  the  value  of  r0  was  significantly  degraded  by 
failing  to  account  for  the  temporal  and  spatial  interaction  of  all  surfaces.  When  only 
accounting  for  a  single  surface  per  detector,  we  found  that  likelihood  was  maximized 
for  a  low  estimate  of  r0.  When  working  with  experimentally  collected  data,  the 
ability  to  accurately  estimate  r0  was  again  degraded.  For  both  target  configurations, 
the  maximization  of  likelihood  with  respect  to  occurred  at  a  value  higher  than 
previously  estimated  or  measured  as  demonstrated  in  Figures  5.19  and  5.20. 


(a)  Adjusted  Likelihood  vs  r0 


(b)  Zoom  View  of  Adjusted  Likelihood 


Figure  5.19:  Target  configuration  #1  single  surface  estimation,  (a)  Adjusted  likelihood 
for  a  range  of  ro  values  from  0.0001m  to  0.01m.  (b)  Maximization  of  likelihood  occurs 
at  0.0016m  where  the  measured  value  was  0.0012m. 


5.5  Chapter  Summary 


Enhancements  in  the  capability  gained  through  3-D  imaging  are  significant.  Cur¬ 
rent  sensors  may  have  limited  military  utility  due  to  the  maturity  of  the  technology. 
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rQ  (meters) 

(a)  Adjusted  Likelihood  vs  r0 


(b)  Zoom  View  of  Adjusted  Likelihood 


Figure  5.20:  Target  configuration  #2  single  surface  estimation,  (a)  Adjusted  likelihood 
for  a  range  of  ro  values  from  0.0001m  to  0.01m.  (b)  Maximization  of  likelihood  occurs 
at  0.0024m  where  the  measured  value  was  0.0018m. 


However,  it  is  expected  that  this  technology  may  eventually  fully  coexist  with  pas¬ 
sive  type  sensors  due  to  a  myriad  of  associated  advantages.  The  algorithm  developed 
for  this  research  demonstrates  that  parameterized  blind  deconvolution  is  an  over¬ 
determined  problem  for  range  diverse  scenes,  with  a  direct  solution  that  maximizes 
likelihood.  When  coupled  with  the  MSID  algorithm,  we  can  simultaneously  discrim¬ 
inate  the  range  to  multiple  surfaces  per  pixel,  while  also  improving  spatial  resolution 
and  temporal  accuracy.  This  algorithm  will  further  enhance  the  ability  to  detect  tar¬ 
gets  behind  an  obscuration.  Improvements  such  as  these  are  critical  for  the  migration 
of  this  technology  to  the  next  generation  of  imaging  sensors. 

This  algorithm  is  novel  in  its  approach  to  the  parameterized  blind  deconvolution 
problem  in  that  it  uses  the  added  information  available  with  3-D  imaging.  Given  the 
following  conditions: 

•  Exposure  time  is  short  enough  that  the  atmosphere  can  be  considered  static 

•  Range  diversity  exists  in  the  targeted  scene 

•  Relationship  in  (5.3)  is  satisfied 
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the  problem  appears  to  have  a  direct  solution.  All  work  to  this  point  has  been  con¬ 
ducted  with  the  assumption  that  the  images  are  properly  sampled.  This  requirement 
was  a  driving  factor  for  using  such  high  levels  of  turbulence  in  the  experiments.  It 
is  possible  that  similar  techniques  could  be  used  for  images  that  are  not  properly 
sampled  or  slightly  under  sampled  and  this  possibility  poses  an  interesting  topic  for 
future  experimentation.  Finally,  while  this  technique  was  only  demonstrated  with  the 
short  exposure  OTF  parameterized  by  ro,  it  is  likely  that  the  technique  would  work 
for  other  blurring  functions  that  can  be  reduced  to  functions  of  a  few  parameters  such 
as  focus  error. 
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VI.  Conclusion 


Contributions  from  this  research  enhance  the  ability  to  accurately  extract  pulse 
information  with  a  three  dimensional  FLASH  LAser  Detection  and  Ranging  (3-D 
FLASH  LADAR)  sensor  where  multiple  returns  per  image  pixel  are  possible.  The 
algorithms  and  techniques  developed  significantly  enhance  ranging  accuracy  and  tar¬ 
get  modeling  where  the  received  image  is  degraded  by  the  effects  of  imaging  through 
a  turbid  medium.  Additionally,  the  algorithms  were  designed  with  computational 
efficiency  in  mind  such  that  they  could  be  employed  in  a  tactical  environment  where 
processing  time  and  power  may  be  limited. 

Active  illumination  sensors  such  as  3-D  FLASH  LADAR  have  recently  garnered 
a  significant  amount  of  interest  for  defense  and  civilian  applications.  Due  to  the 
employment  of  active  flood  illumination,  3-D  FLASH  LADAR  sensors  can  gather 
ranging  information  for  every  point  within  a  targeted  scene  nearly  simultaneously. 
Furthermore,  depending  on  the  detection  methodology  employed,  the  possibility  of 
imaging  through  obscurations  becomes  a  reality.  Potential  benefits  from  this  technol¬ 
ogy  are  currently  constrained  by  manufacturing  challenges  as  well  as  the  difficulty  in 
extracting  the  parameters  of  interest  out  of  the  received  data.  This  research  focuses 
on  the  latter  constraint  of  parameter  extraction  and  enhancement.  As  an  aside,  the 
research  concludes  that  unique  solutions  exist  to  some  problems  that  were  previously 
deemed  mathematically  ill-posed  for  traditional  optical  sensors. 

This  chapter  provides  a  summary  of  each  of  the  primary  chapters  in  this  disser¬ 
tation,  reviews  the  key  contributions  to  this  research  discipline  and  offers  numerous 
ideas  for  future  research  that  could  yield  even  further  enhancement. 
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6.1  Chapter  Summaries 


Chapter  II  reviewed  key  theoretical  concepts  that  were  employed  throughout  this 
research  effort.  The  focus  of  this  review  was  on  3-D  FLASH  LADAR  operation, 
statistical  modeling  of  noise  in  collected  images,  estimation  techniques  with  itera¬ 
tive  algorithms  and  the  effects  of  atmospheric  turbulence  on  collected  imagery.  An 
overview  of  the  sensors  employed  for  collection  of  experimental  data  was  provided.  Fi¬ 
nally,  previous  research  on  blind  deconvolution,  multi-surface  ranging  and  3-D  FLASH 
LADAR  image  enhancement  was  summarized  and  compared  with  this  work. 

Chapter  III  revisited  parameterized  blind  deconvolution  techniques  that  were  de¬ 
veloped  by  MacDonald  and  MacManus  for  two  dimensional  (2-D)  imagery  [41],  [45]. 
MacDonald  employed  a  noise  based  stopping  criteria  in  conjunction  with  an  a  pri¬ 
ori  assumption  for  the  distribution  of  the  atmospheric  seeing  parameter,  r0.  This 
allowed  him  to  simultaneously  find  a  maximum  likelihood  solution  for  the  parame¬ 
terized  Optical  Transfer  Function  (OTF)  and  the  deblurred  image.  MacManus  later 
demonstrated  that  the  a  priori  assumption  was  not  required.  Rather,  employment 
of  the  noise  based  stopping  criteria  was  sufficient  to  solve  the  parameterized  blind 
deconvolution  problem.  Numerous  experiments  with  fully  illuminated  and  partially 
illuminated  scenes  were  conducted  to  validate  the  results  obtained  through  simu¬ 
lation.  Additionally,  the  algorithm  was  employed  against  scenes  illuminated  with 
both  incoherent  and  partially  coherent  light.  Finally,  this  technique  was  applied  to 
simulated  3-D  FLASH  LADAR  images  to  demonstrate  its  effectiveness. 

Chapter  IV  details  the  development  of  the  novel  multiple  surface  ranging  algo¬ 
rithm  that  accounts  for  the  effects  of  diffraction.  This  algorithm  takes  the  somewhat 
common  approach  of  applying  a  Gaussian  mixture  model  to  the  received  3-D  FLASH 
LADAR  data.  However,  the  pulse  model  is  expanded  to  simultaneously  incorporate 
the  effects  of  diffraction  on  the  received  pulse.  Solutions  for  pulse  amplitude,  range 
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and  width  are  developed  for  each  return  in  the  received  data.  Additionally  a  solution 
for  the  detector  bias  is  developed  for  situations  where  it  is  impractical  to  measure 
separately.  A  novel  approach  to  determining  the  optimal  mixture  or  number  of  re¬ 
turns  per  pixel  is  presented.  Finally,  a  constraint  on  the  pulse  amplitude  estimates 
is  derived  to  prevent  improper  convergence  of  the  algorithm.  Simulated  3-D  FLASH 
LADAR  data  was  used  to  compare  the  performance  of  the  various  ranging  techniques. 
Traditional  Gaussian  mixture  modeling  without  considering  the  effects  of  diffraction, 
frame-by-frame  2-D  deconvolution  techniques  and  the  constrained  /  non-constrained 
multi-surface  ranging  algorithms  were  each  employed  on  the  simulated  data. 

Chapter  V  provides  the  capstone  for  the  research  effort.  In  this  chapter  the  set 
of  equations  are  developed  that  demonstrate  that  parameterized  blind  deconvolution 
is  often  an  overdetermined  problem  when  working  with  range  diverse  3-D  FLASH 
LADAR  data.  A  maximum  likelihood  solution  is  employed  to  determine  the  opti¬ 
mal  level  of  atmospheric  seeing  to  deblur  the  collected  images.  The  technique  is 
first  employed  on  simulated  data  where  the  target  geometry,  orientation  and  atmo¬ 
spheric  seeing  are  precisely  controlled.  An  actual  3-D  FLASH  LADAR  sensor  is  then 
employed  to  validate  the  results  from  simulation. 

6.2  Summary  of  Key  Contributions 

6.2.1  Multi-surface  Ranging  Algorithm. 

Perhaps  the  most  significant  contribution  of  this  work  was  the  derivation  of  a 
multi-surface  ranging  algorithm  that  incorporates  the  effects  of  diffraction.  This 
algorithm  simultaneously  estimates  pulse  parameters  from  the  received  data  for  mul¬ 
tiple  surfaces  per  image  pixel.  Through  a  survey  of  relevant  literature,  no  other 
multi-surface  ranging  algorithm  that  removes  the  effects  of  diffraction  could  be  iden¬ 
tified.  Depending  on  the  application  of  the  3-D  FLASH  LADAR  sensor,  the  effects 
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of  diffraction  can  have  significant  negative  impact  on  the  received  images. 

Multiple  enhancements  were  observed  in  the  processed  images  with  this  algorithm. 
First,  the  accuracy  of  this  multi-surface  estimation  is  improved  through  the  reduc¬ 
tion  of  diffraction  effects  caused  by  imaging  through  a  turbulent  atmosphere.  With 
2-D  images,  diffraction  effects  manifest  as  a  spatial  blur.  However,  in  3-D  imaging, 
diffraction  effects  simultaneously  manifest  as  a  spatial  and  temporal  blur.  The  tem¬ 
poral  blurring  can  cause  the  realization  of  false  surface  returns  in  the  received  data. 
The  algorithm  developed  for  this  research  significantly  reduces  temporal  blurring. 
Traditional  multi-surface  ranging  algorithms  would  treat  the  surfaces  that  arise  from 
temporal  blurring  as  additional  returns,  while  the  MSID  algorithm  attempts  to  re¬ 
move  them.  Second,  the  employment  strategy  for  the  algorithm  discussed  within  this 
dissertation  addresses  common  challenges  associated  with  Expectation  Maximization 
(EM)  for  Gaussian  mixture  modeling  such  as  parameter  initialization  and  more  im¬ 
portantly  determining  the  number  of  components  that  exist  in  the  mixture.  This 
research  employs  a  Cap  and  Refine  (CaR)  strategy,  where  the  presence  of  false  re¬ 
turns  are  initially  accounted  for  but  then  eliminated  through  iteration  of  the  EM 
algorithm.  Third,  the  algorithm  also  improves  the  spatial  representation  of  the  image 
that  is  commonly  associated  with  2-D  blind  deconvolution  techniques. 

6.2.2  Amplitude  Constraint  to  Prevent  Early  Convergence  of  Vari¬ 
ance. 

Convergence  of  variance  is  a  useful  stopping  criteria  for  algorithms  such  as  the 
one  developed  in  this  dissertation.  However,  through  the  course  of  testing  this  algo¬ 
rithm,  it  was  uncovered  that  the  algorithm  commonly  terminates  before  the  optimal 
solution  is  found.  This  is  due  to  the  simultaneous  estimation  of  pulse  amplitude  and 
bias  parameters.  A  constraint  was  derived  through  the  use  of  a  Lagrange  multiplier 
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that  significantly  reduced  the  chances  for  improper  convergence.  By  applying  a  con¬ 
straint  on  the  amplitude  that  is  based  on  the  deblurred  2-D  image  representation,  we 
minimize  the  chances  of  hitting  the  stopping  criteria  for  the  iterative  algorithm  before 
the  optimal  solution  is  found.  Of  notable  importance,  this  amplitude  constraint  has 
application  even  in  cases  where  the  image  is  not  impacted  by  the  effects  of  diffraction. 

6.2.3  Maximum  Likelihood  Solution  to  Parameterized  Blind  Decon¬ 
volution  Problem. 

Through  consideration  of  range  diversity  available  with  3-D  FLASH  LADAR,  the 
problem  of  parameterized  blind  deconvolution  often  becomes  over-determined.  Using 
the  multi-surface  ranging  algorithm,  this  research  shows  that  likelihood  is  maximized 
for  the  correct  value  of  atmospheric  seeing.  Due  to  the  ill-posed  nature  of  the  problem 
with  2-D  images,  similar  attempts  at  a  solution  require  additional  assumptions  and/or 
a  priori  information  to  allow  for  a  likelihood  maximization  approach.  Of  additional 
importance,  the  capability  to  produce  a  solution  to  this  previously  ill-posed  problem 
by  taking  advantage  of  range  diversity  highlights  potential  opportunities  for  other 
challenging  problems  commonly  encountered  with  2-D  imaging. 

6.2.4  Employment  of  CoV  Technique  on  3-D  FLASH  LADAR  Images. 

The  results  presented  in  this  dissertation  demonstrate  that  the  Convergence  of 
Variance  (CoV)  technique  originally  developed  by  MacManus  has  application  in  de¬ 
termining  the  parameterized  blurring  function  in  3-D  FLASH  LADAR  imaging.  How¬ 
ever,  this  research  also  points  to  the  relationship  between  accuracy  of  the  r0  estimates 
and  the  Signal  to  Noise  Ratio  (SNR)  in  the  collected  images.  Due  to  the  manner  in 
which  3-D  FLASH  LADAR  images  are  recorded,  the  user  may  be  required  to  derive 
an  optimal  method  for  selection  of  a  2-D  frame  on  which  the  CoV  technique  will  be 
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performed. 


6.3  Future  Research 

During  the  execution  of  this  research  program,  numerous  topics  were  uncovered 
which  may  yield  further  enhancement  in  the  field  of  imaging  with  3-D  FLASH  LADAR 
sensors. 

6.3.1  2-D  and  3-D  Data  Fusion. 

The  algorithms  developed  in  this  research  effort  and  the  results  they  produced 
relied  on  properly  sampled  imaging  data.  However,  the  current  state  of  sensor  tech¬ 
nology  presents  numerous  challenges  to  the  collection  of  properly  sampled  data.  This 
research  effort  obtained  properly  sampled  data  by  using  an  extremely  turbulent  at¬ 
mosphere  to  act  as  a  low-pass  filter.  Alternatively,  given  the  current  state  of  the 
technology,  properly  sampled  data  could  only  be  achieved  with  optical  configurations 
with  high  focal  ratios.  Through  the  fusion  of  properly  sampled  2-D  imagery  and  under 
sampled  3-D  imagery,  it  may  be  possible  to  reduce  the  sampling  requirement  for  3-D 
sensors  while  yielding  similar  performance  gains  summarized  in  this  work.  We  found 
that  the  multi-surface  ranging  algorithm  had  some  capability  to  deal  with  slightly 
undersampled  data.  However,  the  extent  to  which  the  data  can  be  undersampled 
with  this  technique  also  remains  to  be  proven  as  a  future  research  topic. 

6.3.2  Applicability  to  Other  Parameterized  Blurring  Functions. 

For  purposes  of  this  research,  we  substituted  the  average  short  exposure  transfer 
function,  Hse,  as  our  model  for  the  atmospheric  Optical  Transfer  Function  (OTF), 
Hatm.  However,  the  techniques  developed  should  be  directly  applicable  to  other  cases 
where  the  blurring  function  can  be  reduced  to  relatively  few  unknowns.  One  inter- 
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esting  topic  would  be  to  extend  this  research  to  remove  the  effects  of  focus  blur.  This 
particular  extension  may  be  useful  for  deblurring  3-D  FLASH  LADAR  images  with 
long  range  gates  at  relatively  short  ranges.  Additionally,  a  potential  extension  for  this 
future  research  would  be  the  relationship  between  the  seeing  parameter  and  focus  in¬ 
teraction  for  varying  levels  of  focus  error  and  atmospheric  turbulence.  In  Chapter  III 
we  noticed  an  interaction  between  focus  and  atmospheric  blur  which  resulted  in  low 
estimates  for  atmospheric  seeing.  It  would  be  interesting  to  see  if  the  two  could  be 
separated. 

6.3.3  Direct  Solution  for  Parameterized  OTF. 

This  research  employed  a  maximization  of  likelihood  approach  to  find  the  optimal 
parameterized  OTF.  However,  it  may  be  possible  to  directly  solve  for  the  parameter¬ 
ized  OTF  in  conjunction  with  the  deblurred  pulse  model.  An  attempt  was  made  to 
approximate  the  Point  Spread  Function  (PSF)  as  a  Gaussian  and  directly  solve  for 
the  width  parameter  in  an  iterative  fashion.  However,  this  technique  did  not  yield 
accurate  results. 

6.3.4  Alternative  Stopping  Criteria. 

The  iterative  algorithms  presented  in  this  research  rely  on  the  CoV  technique  to 
establish  the  stopping  criteria.  However,  due  to  the  Gaussian  approximation  for  the 
returned  pulse,  this  technique  may  be  problematic  in  non-uniform  data.  An  explo¬ 
ration  /  comparison  of  other  stopping  criteria  is  warranted  such  as  those  presented 
by  Arioli  [3]. 
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